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We focus on mesoscopic dislocation patterning via a continuum dislocation dynamics theory 
(CDD) in three dimensions (3D). We study three distinct physically motivated dynamics which 
consistently lead to fractal formation in 3D with rather similar morphologies, and therefore we sug- 
gest that this is a general feature of the 3D collective behavior of geometrically necessary dislocation 
(GND) ensembles. The striking self-similar features are measured in terms of correlation functions 
of physical observables, such as the GND density, the plastic distortion, and the crystalline orien- 
tation. Remarkably, all these correlation functions exhibit spatial power-law behaviors, sharing a 
single underlying universal critical exponent for each type of dynamics. 

PACS numbers: 61.72.Bb, 61.72.Lk, 05.45.Df, 05.45.Pq 

I. INTRODUCTION 

Dislocations in plastically deformed crystals, driven by their long-range interactions, collectively evolve into complex 
heterogeneous structures where dislocation-rich cell walls or boundaries surround dislocation-depleted cell interiors. 
These have been observed both in single crystal^^^^^ and polycrystalP^ using transmission electron microscopy (TEM) . 
The mesoscopic cellular structures have been recognized as scale-free patterns through fractal analysis of TEM mi- 
crograph^^HH xhe complex collective behavior of dislocations has been a challenge for understanding the underlying 
physical mechanisms responsible for the development of emergent dislocation morphologies. 

Complex dislocation microstructures, as an emergent mesoscale phenomenon, have been previously modeled using 
various theoretical and numerical approacheJ^. Discrete dislocation dynamics (DDD) models have provided insights 
into the dislocation pattern formations: parallel edge dislocations in a two-dimensional system evolve into 'matrix 
structures' during single slip^*^, and 'fractal and cell structures' during multiple slip^^ random dislocations in a three- 
dimensional syst em se lf-organize themselves into microstructures through junction formation, cross-slip, and short- 
range interactionJl2lIII However, DDD simulations are limited by the computational challenges on the relevant scales 
of length and strain. Beyond these micro-scale descriptions, CDD has also been used to study complex dislocation 
structures. Simplified reaction-diffusion models have described persistent slip bands^^, dislocation cellular structures 
during multiple slipf^Sl, and dislocation vein structurefP^ . Stochasticity in CDD models^ loH! q,. jj^ ^j^g splittings 
and rotations of the macroscopic cellJ^^t^Il have been suggested as an explanation for the formation of organized 
dislocation structures. The source of the noise in these stochastic theories is derived from cither extrinsic disorder or 
short-length-scale fluctuations. 

In a recent manuscriplpS, we analyzed the behavior of a grossly simplified continuum dislocation model for plastic- 
itjl22H26|_ physicist's 'spherical cow' approximation designed to explore the minimal ingredients necessary to explain 
key features of the dynamics of deformation. Our simplified model ignores many features known to be important for 
cell boundary morphology and evolution, including slip systems and crystalline anisotropy, dislocation nucleation, lock 
formation and entanglement, line tension, geometrically unnecessary forest dislocations, etc. H owever, our model does 
encompass a realistic order parameter field (the Nye-Kroner dislocation density tensoi'22l2l! embodying the GNDs) , 
which allows detailed comparisons of local rotations and deformations, stress, and strain. It is not a realistic model 
of a real material, but it is a model material with a physically sensible evolution law. Given these simplifications, 
our model exhibited a surprisingly realistic evolution of cellular structures. We analyzed these structures in two- 
dimensional simulations (full three-dimensional rotations and deformations, but uniform along the z-axis) using both 
the fractal box counting method^ * and the single-length-scale scaling methodJ^SMH used in previous theoretical anal- 
yses of experimental data. Our model qualitatively reproduced the self-similar, fractal patterns found in the former, 
and the scaling behavior of the cell sizes and misorientations under strain found in the latter (power-law refinement 
of the cell sizes, power-law increases in misorientations, and scaling collapses of the distributions). 

There are many features of real materials which are not explained by our model. We do not observe distinctions 
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FIG. 1: Experimental and simulated dislocation cellular structures. In (a), a typical TEM micrograph at a micron 
scale is talten from a Cu single crystal after [100] tensile deformation to a stress of 76.5 MPar* in (b), a simulated GND density 
plot is shown. Note the striking morphological similarity between theory and experiment. 



between 'geometrically necessary' and 'incidental' boundaries, which appear experimentally to scale in different ways. 
The fractal scaling observed in our model may well be cut off or modified by entanglement, slip-system physics, 
quantization of Burgers vector'^ or anisotropy - we cannot predict that real materials should have fractal cellular 
structures; we only observe that our model material does so naturally. Our spherically symmetric model obviously 
cannot reproduce the dependence of morphological evolution on the axis of app lied strain (and hence the number of 
activated slip systems); indeed, t he frac tal patterns observed in some experiments^ could be associated with the high- 
symmetry geometry they studiecP^I^. While many realistic features of materials that we ignore may be important 
for cell-structure formation and evolution, our model gives clear evidence that these features are not essential to the 
formation of cellular structures when crystals undergo plastic deformation. 

In this longer manuscript, we provide an in-depth analysis of three plasticity models. We show how they (and more 
traditional models) can be derived from the structures of the broken symmetries and order parameters. We extend 
our simulations to 3D, where the behavior is qualitatively similar with a few important changes. Here we focus our 
attention on relaxation (rather than strain), and on correlation functions (rather than fractal box counting or cell 
sizes and misorientations) . 

Studying simplified 'spherical cow' models such as ours is justified if they capture some key phenomenon, providing 
a perspective or explanation for the emergent behavior. Under some circumstances, these simplified models can 
capture the long-wavelength behavior precisely - the model is said to be in the same universality class as the observed 
behavior^ (Chapter 12). The Ising model for magnetism, two- fluid criticality, and order-disorder transitions; self- 
organized critical models for magnetic Barkhausen noise^^ '^^ and dislocation avalancheS^ all exhibit the same type 
of emergent scale-invariant behavior as observed in some experimental cellular structureJ^. For all of these systems, 
'spherical cow' models provide quantitative experimental predictions of all phenomena on long length and time scales, 
up to overall amplitudes, relevant perturbations, and corrections to scaling. Other experimental cellular structureP^l 
have been interpreted in terms of scaling functions with a characteristic scale, analogous to those seen in crystalline 
grain growth. Crystalline grain growth also has a 'universal' description, albeit one which depends upon the entire 
anisotropic interfacial energy and mobilit}'^ (and not just temperature and field) .1221 We are cautiously optimistic that 
a model like ours (but with metastability and crystalline) could indeed describe the emergent complex dislocation 
structures and dynamics in real materials. Indeed, recent work on dislocation avalanches suggests that even the yield 
stress may be a universal critical poinlP^l. 

Despite universality, we must justify and explain the form of the CDD model we study. In Sec. |lT] we take the 
continuum, 'hydrodynamic' limit approach, traditionally originating with Landau in the study of systems near thermal 
equilibrium (clearly not true of deformed metals!). All degrees of freedom are assumed slaves to the order parameter, 
which is systematically constructed from conserved quantities and broken symmetrieS^^HMI _ t^jg jg the fundamental 
tool used in the physics community to derive the diffusion equation, the Navier-Stokes equation, and continuum 
equations for superconductors, superfluids, liquid crystals, etc. Ref . H51 have utilized this general approach to generate 
CDD theories, and in Sec. |ll] we explain how our approach differs from theirs. 

In Sec. |III| we explore the validity of several approximations in our model, starting in the engineering language of 
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state variables. Here local equilibration is not presumed; the state of the system depends in some arbitrarily complex 
way on the history. Conserved quantities and broken symmetries can be supplemented by internal state variables - 
statistically stored dislocations (SSDs), yield surfaces, void fractions, etc., whose evolution laws are judged according 
to their success in matching experimental observations. (Eddy viscosity theories of turbulence are particular successful 
examples of this framework.) The 'single-velocity' models we use were originally developed by Acharya et al?"^ 
and we discuss their microscopic derivatiorP^ and the correction term resulting from coarse-graining and multiple 
microscopic velocitie^2l!_ This term is usually modeled by the effects of SSDs using crystal plasticity models. We 
analyze experiments to suggest that ignoring SSDs may be justified on the length-scales needed in our modeling. 
However, we acknowledge the near certainty that Acharya's will be important - the true coarse-grained evolution 
laws will incorporate multiple velocities. Our model should be viewed as a physically sensible model material, not a 
rigorous continuum limit of a real material. 

In this manuscript, we study fractal cell structures that form upon relaxation from randomly deformed initial 
conditions (Sec. IV B ). One might be concerned that relaxation of a randomly imposed high-stress dislocation structure 
(an instantaneous hammer blow) could yield qualitatively different behavior from realistic deformations, where the 
dislocation structures evolve continuously as the deformation is imposed. In Sec. |IVB| we note that this alternative 
'slow hammering' gives qualitatively the same fractal dislocation patterns. Also, the re sultin g cellular structures are 
qualitatively very similar to those we observe under subsequent uniform external straiiP^'^, except that the relaxed 
structures are statistically isotropic. We also find that cellular structures form immediately at small deformations. 
Cellular structures in real materials emerge only after significant deformation; presumably this feature is missing in 
our model because our model has no impediment to cross-slip or multiple slip, and no entanglement of dislocations. 
This initial relaxation should not be viewed as annealing or dislocation creep. A proper description of annealing 
must include dislocation line tension effects, since the driving force for annealing is the reduction in total dislocation 
density - our dislocations annihilate when their Nye Burgers vector density cancels under evolution, not because of 
the dislocation core energies. Creep involves dislocation climb, which (for two of our three models) is forbidden. 

We focus here on correlation functions, rather than the methods used in previous analyses of experiments. Correla- 
tion functions have a long, dignified history in the study of systems exhibiting emergent scale invariance - materials at 
continuous thermodynamic phase transitions'*-^, fully developed turbulencc*^^^ and crackling noise and self-organized 
criticalitjEH. We study not only numerical simulations of these correlations, but provide also extensive analysis of the 
relations between the correlation functions for different physical quantities and their (possibly universal) power-law 
exponents. The decomposition of the system into cells (needed for the cell-size and misorientation distribution anal- 
ygg^29H32lj demands the introduction of an artificial cutoff misorientation angle, and demands either laborious human 
work or rather sophisticated numerical algorithm^^. These sections of the current manuscript may be viewed both 
as a full characterization of the behavior of our simple model, and as an illustration of how one can use correlation 
functions to analyze the complex morphologies in more realistic models and in experiments providing 2D or 3D real- 
space data. We believe that analyses that explicitly decompose structures into cells remain important for systems 
with single changing length-scale: grain boundary coarsening should be studied both with correlation functions and 
with explicit studies of grain shape and geometry evolution, and the same should apply to cell-structure models and 
experiments that are not fractal. But our model, without such an intermediate length-scale, is best analyzed using 
correlation functions. 

Our earlier worlP^l focused on 2D. How different are our predictions in 3D? In this paper, we explore three different 
CDDs that display similar dislocation fractal formation in 3D and confirm analytically that correlation functions of 
the GND density, the plastic distortion, and the crystalline orientation, all share a single underlying critical exponent, 
up to exponent relations, dependent only on the type of dynamics. Unlike our 2D simulations, where forbidding climb 
led to rather distinct critical exponents, all three dynamics in 3D share quite similar scaling behaviors. 

We begin our discussion in Sec. |II A| by defining the various dislocation, distortion, and orientation fields. In Sec. |II B[ 
we derive standard local dynamical evolution laws using traditional condensed matter approaches, starting from both 
the non-conserved plastic distortion and the conserved GND densities as order parameters. Here, we also explain 
why these resulting dynamical laws are inappropriate at the mesoscale. In Sec. |II C[ we show how to extend this 
approach by defining appropriate constitutive laws for the dislocation flow velocit y to build no vel dynamics^^. There 
are three different dynamics we study: i) isotropic climb-and-glide dynamics (CGD)2 3 | 24 | 26 | 53 | 5 3 jj) isotropic glide-only 
dynamics, where we define the part of the local dislocation density that participates in the local mobile dislocation 
population, keeping the local volume conserved at all times (GOD-MDP)^^, iii) isotropic glide-only dynamics, where 
glide is enforced by a local vacancy pressure due to a co-existing background of vacancies that have an infinite energy 
cost (G0D-LVP)2S1. All three types of dynamics present physically valid alternative approaches for deriving a coarse- 



grained continuum model for GNDs. In Sec. Ill we explore the effects of coarse-graining, explain our rationale 
for ignoring SSDs at the mesoscale, and discuss the single-velocity approximation we use. In Sec. |IV[ we discuss the 
details of numerical simulations in both two and three dimensions, and characterize the self-organized critical complex 
patterns in terms of correlation functions of the order parameter fields. In Sec. [Vj we provide a scaling theory, derive 
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relations among the critical exponents of these related correlation functions, study the correlation function as a scaling 
function of coarse-graining length scale, and conclude in Sec. |VI| 

In addition, we provide extensive details of our study in Appendices. In [Xj we collect useful formulas from the 
literature relating different physical quantities within traditional plasticity, while in[B]we show how functional deriva- 
tives and the dissipation rate can be calculated using this formalism, leading to our proof that our CDDs are strictly 
dissipative (lowering the appropriate free energy with time). In[C| we show the flexibility of our CDDs by extending 
our dynamics: In particular, we show how to add vacancy diffusion in the structure of CDD, and also, how external 
disorder can be in principle incorporated (to be explored in future work). In|D] we elaborate on numerical details - we 
demonstrate the statistical convergence of our simulation method and also we explain how we construct the Gaussian 
random initial conditions. Finally, in [Ej we discuss the scaling properties of several correlation functions in real 
and Fourier spaces, including the strain- history-dependent plastic deformation and distortion fields, the stress-stress 
correlation functions, the elastic energy density spectrum, and the stressful part of GND density. 



II. CONTINUUM MODELS 
A. Order parameter fields 

1. Conserved order parameter field 

A dislocation is the topological defect of a crystal lattice. In a continuum theory, it can be described by a coarse- 
grained variable, the GND density,!^ (also called the net dislocation density or the Nye-Kroner dislocation density), 
which can be defined by the GND density tensor 

p(x) = ^(t".n)n0b"5(x-r), (1) 



PkrrM)^Y.^tb^^^^-n, (2) 

a 

measuring the sum of the net flux of dislocations a located at ^, tangent to t, with Burgers vector b, in the neigh- 
borhood of x, through an infinitesimal plane with the normal direction along n, seen in Fig. [2j In the continuum, the 
discrete sum of line singularities in Eqs. ([l]) and ^ is smeared into a continuous (nine-component) field, just as the 
continuum density of a liquid is at root a sum of point contributions from atomic nuclei. 

Since the normal unit pseudo- vector fi is equivalent to an antisymmetric unit bivector Eij — Sijkfik, we can 
reformulate the GND density as a three-index tensor 

e(x) = V(t".n)^^b"<5(x~r), (3) 



SO 



Qiji 



,(x)=5](t".n)^„-0(x-r), (4) 



measuring the same sum of the net fiux of dislocations in the neighborhood of x. t hrough the infinitesimal plane 



indicated by the unit bivector E. This three-index variant will be useful in Sec. II C 2 where we adapt the equations 
of Refs.i2i]and[lS]to forbid dislocation climb (GOD-MDP). 

According to the definition of E, we can find the relation between p and g 

Qi]m{^) = ^{i'lfli)eijknkKn^{'^ - = SijkPkmi^)- (5) 

a 

It should be noted here that dislocations cannot terminate within the crystal, implying that 

9iPy(x) = 0, (6) 

or 



eijkdkQijiix) = 0. 



(7) 




FIG. 2: (Color online) Representation of the crystalline line defect — dislocation. Each curved line represents a 
dislocation line with the tangent direction t, and the Burgers vector b which characterizes the magnitude and direction of the 
distortion to the lattice. The two-index GND density pk m^"^'^^ ' (Eqs.jljand [2]) is the net flux of the Burgers vector density b 
alonge'-™' through an infinitesimal piece of a plane with normal direction n along e'*^'. The three-index version Qijm (Eqs. [s] 
and k| is the flux density through the plane along the axes e'*' and e^-*', with the unit bivector E — e*-'' A e*--''. 



Within plasticity theories, the gradient of the total displacement field u represent s the compatible total distortion 
j^g y28 | 55 | p,, _ Q.y:.^ which is the sum of the elastic and the plastic distortion field^^HHII^ fj = -\- p^. Due to the 
presence of dislocation lines, both /?p and /3° are incompatible, characterized by the GND density p 



The elastic distortion field is the sum of its symmetric strain and antisymmetric rotation fields, 

= e° + w^ 



(8) 
(9) 



(10) 



where we assume linear elasticity, ignoring the 'geometric nonlinearity' in these tensors. Substituting the sum of two 
tensor fields into the incompatibility relation Eq. ([s]) gives 

Pi j = £iki dkUJij +Sikidkeij. (11) 
The elastic rotation tensor can be rewritten as an axial vector, the crystalline orientation vector A 

1 



Afc 



(12) 



or 



(13) 
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Thus we can substitute Eq. (131 into Eq. (11) 



Pij = [kjdkt^k - djA,) + emdkf^ly (14) 

For a system without residual elastic stress, the GND density thus depends only on the varying crystalline orienta- 
tion^el. ' 

Dynamically, the time evolution law of the GND density emerges from the conservation of the Burgers vectoi'^^'^ 

d 



or 



dt 



Qijk 



%^mpq^pJqk — Qijpq^pJqkf 



where J represents the Burgers vector flux, and the symbol Qijpq indicates 

^ijm^mpq — ^ip^jq ^iq^jp- 



2. N on- conserved order parameter field 



(16) 



The natural physicist's order parameter field g, characterizing the incompatibility, can be written in terms of the 
plastic distortion field (3^ 



Qijk 



iPvak 



-gtjisdiP^k- 



(17) 



In the linear approximation, the alternative order parameter field (3^ fully specifies the local deformation u of the 
material, the elastic distortion (3°, the internal long-range stress field cr'"* and the crystalline orientation (the Rodrigues 
vector A giving the axis and angle of rotation) , as summarized in [A 

It is natural, given Eq. ^ and Eq. to use the flux J of the Burgers vector density to define the dynamics of 
the plastic distortion tensor ^nlMZiMt 



dt 



Ji' 



(18) 



As noted by Ref. [Ml Eq. ^ and Eq. (15 1 equate a curl of /3p to a curl of J, so an arbitrary divergence may be 
added to Eq. (18): the evolution of the plastic distortion /?p is not determined by the evolution of the GND density. 
Ref. ,54, resolves this ambiguity using a Stokes-Helmholtz decomposition of (3^. In our notation, /Jp — /3P'^ 



p,H 



The 'intrinsic' plastic distortion /3P'^ is divergence-free {di(3f^ = 0, i.e., hPf: = 0), and determined by the GND 

density p. The 'history- dependent!^ 
decomposition explicitly, as 



P'^ is curl-free {eujdePfj — 0, eujkiPf: = 0). In Fourier space, we can do this 



(k) 



iSilm-j^Pmjik) + iki1pj{]&) 



^l^''(k)+/3r(k) 



P,H/ 



(19) 



This decomposition will become important to us in Sec. IV C 3 where the correlation functions of /3P'^ and /3P'^ will 
scale differently with distance. 

Ref . [51] treats the evolution of the two components /3P'^ and /3P'^ separately. Because our simulations have periodic 
boundary conditions, the evolution of /3P'^ does not affect the evolution of p. As noted by Ref. [511 in more general 
situations /3P'^ will alter the shape of the body, and hence interact with the boundary conditiongSS. Hence in the 



simulations presented here, we use Eq. (18), with the warning that the plastic deformation fields shown in the figures 
are arbitrary up to an overall divergence. The correlation functions we study of the intrinsic plasti c dis tortion /3P'^ 
are independent of this ambiguity, but the correlation functions of /3P'^ we discuss in the Appendix E 1 will depend 
on this choice. 

In the presence of external loading, we can express the appropriate free energy T as the sum of two terms: the elastic 
interaction energy of GNDs, and the energy of interaction with the applied stress field. The free energy functional is 



int c _cxt p 

2 iJ iJ 



Alternatively, it can be rewritten in Fourier space 



(27r)3 



2 -^^2 J mn V 



(k)/3f, (k)/3P „(~k) + ;?,f(k)/3P,(-k) 



(20) 



(21) 



as discussed in IB 1 
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B. Traditional dissipative continuum dynamics 



There are well known approaches for deriving continuum equations of motion for dissipative systems, which in this 
case produce a traditional von Mises-style theorjEl, useful at longer scales. We begin by reproducing these standard 
equations. 

For the sake of simplicity, we ignore external stress (cry- simplified to <j"f) in the following three subsections. We 
start by using the standard methods applied to the non-conserved order parameter , and then turn to the conserved 
order parameter g. 



1. Dissipative dynamics built from the non-conserved order parameter field 

The plastic distortion /?p is a non-conserved order parameter field, which is utilized by the engineering community 
to study texture evolution and plasticity of mechanically deformed structural materials. The simplest dissipative 
dynamics in terms of /Jp minimizes the free energy by steepest descents 

where F is a positive material-dependent constant. We may rewrite it in Fourier space, giving 

The functional derivative 6J^ /6(3f-{—'k) is the negative of the long-range stress 

, , = -A/.,,n„(k)^P„(k) ^ -^?.,(k). (24) 

This dynamics implies a simplified version of von Mises plasticity 

d 



/3f,(k)=Fa,,(k). (25) 



2. Dissipative dynamics built from the conserved order parameter field g 



We can also derive an equation of motion starting from the GND density g, as was done by Ref. For this 
dissipative dynamics Eq. (16), the simplest expression for J is 



Jqk — —^'ablq^l 



5F 

Sgabk ' 



(26) 



where the material-dependent constant tensor F' must be chosen to guarantee a decrease of the free energy with time. 
The infinitesimal change of J- with respect to the GND density g is 



ST[g] 



^3 ST ^ 

d X og. 



5g, 



ijk 



ijk- 



The free energy dissipation rate is thus ST/5t for Sg ~ ^St, hence 

d 



dt 



F[g] 



at 



3 6T dgijk 



Sgijk dt 



Substituting Eq. (16) into Eq. (28 1 and integrating by parts gives 



dt 



J^[g] = Jd^x (^gijpqdp^^^Jqk- 



(27) 



(28) 



(29) 
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Substituting Eq. ( |26[ ) into Eq. (|29| gives 



d_ 



ST 



ST 

Sgabk 



(30) 



Now, to guarantee that energy never increases, we choose rj^^;^ = Tgabiq, (r is a positive material-dependent constant), 
which yields the rate of change of energy as a negative of a perfect square 



ST_ 

bk 



Using Eqs. (16) and (261, we can write the dynamics in terms of g 

d 



ST 



„ Qi3k - ^gi]pq9ablqdpdl - . 
Ot Ogabk 



(31) 



(32) 



Substituting the functional derivative ST /Sgabk, Eq- (BIO I, derived in B2 into Eq. (32) and comparing to Eq. (16) 
tells us 



d_ 
di 



where 



Jqk = TCF 



qk 



(33) 



(34) 



duplicating the von Mises law (Eq. 25 ) of the previous subsection. The simplest dissipative dynamics of either non- 
conserved or conserved order parameter fields thus turns out to be the traditional linear dynamics, a simplified von 
Mises law. 

The problem with this law for us is that it allows for plastic deformation in the absence of dislocations, i.e., the 
Burgers vector flux can be induced through the elastic loading on the boundaries, even in a defect-free medium. This 
is appropriate on engineering length scales above or around a micron, where SSDs dominate the plasti c deformation. 
(Methods to incorporate their effects into a theory like ours have been provided by Acharya et alW^^ and Varadhan 
et al^ 

By ignoring the SSDs, our theory assumes that there is an intermediate coarse-grain length scale, large compared 
to the distance between dislocations and small compared to the distance where the cancelling of dislocations with 
different Burgers vectors dominates the dynamics, discussed in Sec. |III| We believe this latter length scale is given by 
the distance between cell w alls (a s discussed in Sec. IV B). The cell wall misorientations are geometrically necessary. 
On the one hand, it is knowrPilSllthat neighboring cell walls often have misorientations of alternating signs, so that on 
coarse-grain length scales just above the cell wall separation one would expect explicit treatment of the SSDs would 
be necessary. On the other hand, the density of dislocations in cell walls is high, so that a coarse-grain lengl^ much 
smaller than the interesting structures (and hence where we believe SSDs are unimportant) should be possibl^^^J. (Our 
cell structures are fractal, with no characteristic 'cell size'; this coarse-grain length sets the minimum cutoff scale of 
the fractal, and the grain size or inhomogeneity length will set the maximum sca le.) With this assumption, to treat 
the formation of cellular structures, we turn to theories of the form given in Eq. ( 15 ), defined in terms of dislocation 
currents J that depend directly on the local GND density. 



Our CDD model 



The microscopic motion of a dislocation under external strain depends upon temperature. In general, it moves 
quickly along the glide direction, and slowly (or not at all) along the climb direction where vacancy diffusion must 
carry away the atoms. The glide speed can be limited by phonon drag at higher temperatures, or can accelerate 
to nearly the speed of sound at low temperature^^. It is traditional to assume that the dislocation velocity is 
over-damped, and proportional to the component of the force per unit dislocation length in the glide plane. 

To coarse-grain this microscopies, for reasons described in Sec. |III| , we choose a CDD model whose dislocation 
currents vanish when the GND density vanishes, without considering SSDs. Ref. [26| derived a dislocation current J 
for this case using a closure approximation of the underlying microscopies. Their work reproduced (in the case of 
both glide and climb) an earlier dynamical model proposed by Acharya et alT^'^^, who also inc orporate the effects of 
SSDs. We follow the general approach of Acharya and collaborators^^ '^^'^SiSMOj gg^,^ II C 1 to derive an evolution 
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law for dislocations allowed both to glide and climb, and then modify it to remove climb in Sec. II C 2 We derive a 
second variant of glide-only dynamics in Sec. |II C"3| by coupling climb to vacancies and then taking the limit of infinite 
vacancy energy, which reproduces a model proposed earlier by Ref. 25 

In our CGD and GOD-LVP dynamics (Sections II C 1 and II C 3 below), all dislocations in the infinitesimal volume 
at X are moving with a common velocity i'(x). We discuss the validity of this single- velocity form for the equations 
of motion at length in Sec. |III[ together with a discussion of the coarse-graining and the emergence of SSDs. We view 
our simulations as physically sensible 'model materials' - perhaps not the correct theory for any particular material, 
but a sensible framework to generate theories of plastic deformation and explain generic features common to many 
materials. 



1. Climb-glide dynamics (CGD) 

We start with a model presuming (perhaps unphysically) that vacancy diffusion is so fast that dislocations climb 
and glide with equal mobility. The elastic Peach-Koehler force due to the stress (t(x) on the local GND density is 
given by = amkQumk- We assume that the velocity v cx f^^, giving a local constitutive relation 

Vu OC CFrnkQumk- (35) 

How should we determine the proportionality constant between velocity and force? In experimental systems, this 
is complicated by dislocation entanglement and short-range forces between dislocations. Ignoring these features, the 
velocity of each dislocation should depend only on the stress induced by the other dislocations, not the local density 
of dislocations'^. We can incorporate this in an approximate way by making the proportionality factor in Eq. ( 35 ) 



inversely proportional to the GND density. We measure the latter by summing the square of all components of p, 
hence \q\ = \/ QijkQijk/2 and u„ — j^^amkQumk, where D \s a. positive material-dependent constant. This choice has 
the additional important feature that the evolution of a sharp domain wall whose width is limited by the lattice cutoff 
is unchanged when the lattice cutoff is reduced. 
The flux J of the Burgers vector is thu^^ 

J^,=v^g^,,^^,^,eurakQu.r (36) 

Notice that this dynamics satisfies our criterion that J = when there are no GNDs (i.e. , g = 0). Notice also that 
we do not incorporate the effects of SSDs (Acharya's L^^; we discuss this further in Sec. 



Substituting this flux J (Eq. k56|) into the free energy dissipation rate (Eq. B16) gives 



III 



d^x (JijJ^j = - / d^x < 0. (37) 



Details are given in|B 3| 



2. Clide-only dynamics: mobile dislocation population ( GOD-MDP) 



When the temperature is low enough, dislocation climb is negligible, i.e., dislocations can only move in their glide 
planes. Fundamentally, dislocation glide conserves the total number of atoms, which leads to an unchanged local 
volume. Since the local volume change in time is represented by the trace Ja of the flux of the Burgers vector, 
conservative motion of GNDs demands Ju — 0. Ref. derived the equation of motion for dislocation glide only, by 



removing the trace of J from Eq. (36). However, their dynamics fails to guarantee that the free energy monotonically 



decreases. Here we present an alternative approach. 



We can remove the trace of J by modifying the first equality in Eq. (36), 



J' 



Quij 



SijQukk 



(38) 



where = Quij — \5ijQukk can be viewed as a subset of 'mobile' dislocations moving with velocity v' . 
Substituting the current (Eq. 38 1 into the free energy dissipation rate (Eq. B16l gives 



'dt 



''ij i^u Quij)- 



(39) 
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If we choose the velocity v'^ cx a-ijg'^^j, the appropriate free energy monotonically decreases in time. We thus express 



D J 



Y^Quij'^ijj where D is a positive material-dependent constant, and the prefactor l/\g\ is added for the same 



reasons, as discussed in the second paragraph of Sec. |IIC 1| 
The current J' of the Burgers vector is thus writteiP^' 



(40) 



This natural evolution law becomes much less self-evident when expressed in terms of the traditional two-index version 
p(Eqs.[l]^ 

r, ^ D_f _ _1 1 

^ij — \^inPmnpmj ^mnpinpmj ^^mmpniPnj H" ^^7nmpinpnj 



( _ _i 1 "l^ 

y^knPmnPmk ^mnPknPmk ^^mmPnkPnk H~ ^^mmPknPnk J j ■> 



(41) 



(which is why we introduce the three- index variant g). 

This current J' makes the free energy dissipation rate the negative of a perfect square in Eq. (B18). Details are 
given in |B 3| 



3. Glide- only dynamics: local vacancy-induced pressure (GOD-LVP) 



At high temperature, the fast vacancy diffusion leads to dislocation climb out of the glide direction. As the 
temperature decreases, vacancies are frozen out so that dislocations only slip in the glide planes. In |C H we present a 
dynamical model coupling the vacancy diffusion to our CDD model. Here we consider the limit of frozen-out vacancies 
with infinite energy costs, which leads to another version of glide-only dynamics. 

According to the coupling dynamics Eq. ( C8 ) , wc write down the general form of dislocation current 



J'' 



where p is the local pressure due to vacancies. 



D 

H 



I Quij t 



(42) 



The limit of infinitely costly vacancies (a — oo in C 1 1 leads to the traceless current, J[[ = 0. Solving this equation 
gives a critical local pressure p"^ 



P 



^pqQspqQskk 



Quaa Qubh 

The corresponding current J" of the Burgers vector in this hmit is thus written 

^pq QspqQskk 



„ _ D 



Quaa Qubh 



mn I Qumn Quij i 



(43) 



(44) 



reproducing the glide-only dynamics proposed by Ref. 25 , 



Substituting the current (Eq. 44 1 into the free energy dissipation rate (Eq. B16l gives 



D 



PK fPK 



in, 



d^f, 



PK\ 2-1 



<o. 



(45) 



where ff^ — (JmnQimn and di — gikk- The equality emerges when the force is along the same direction as d. 

Unlike the traditional linear dissipative models, our CDD model, coarse grained from microscopic interactions, 
drives the random plastic distortion to non-trivial stress-free states with dislocation wall singularities, as schematically 
illustrated in Fig. |3] 

Our minimal CDD model, consisting of GNDs evolving under the long-range interaction, provides a framework 

for understanding dislocation morphologies at the mesosc ale. E ventually, it can be extended to include vacancies by 

coupling them to the dislocation current (as discussed in C 1 or extended to include disorder, dislocation pinning, 

and entanglement by adding appropriate interactions to the free energy functional and refining the effective stress 

field (as discussed in C2|. It has already been extended to include SSDs incorporating traditional crystal plasticity 
theorie&25 59 60.^ 
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Random initial state 



lit 




FIG. 3: (Color online) Relaxation of various CDD models. The blue dot represents the initial random plastically-deformed 
state; the red dots indicate the equilibrated stress-free states driven by different dynamics. Curve A: steepest decent dynamics 
leads to the trivial homogeneous equilibrated state, discussed in Sec. 



|II B[ Curve B: our CDD models settle the system into 
non-trivial stress-free states with wall-like singularities of the GND density, discussed in Sec. |II C| 



III. COARSE GRAINING 



The discussion in Sec. [11] uses the language and conceptual framework of the condensed matter physics of systems 
close to equilibrium - the generalized "hydrodynamics" used to derive equations of motion for liquids and gases, liquid 
crystals, superfluids and superconductors, magnetic materials, etc. In these subjects, one takes the broken symmetries 
and conserved quantities, and systematically writes the most general evolution laws allowed by symmetry, presuming 
that these quantities determine the state of the material. In that framework, the Burgers vector flux J of Eqs. ( 15 ) 



and ( 16 ) would normally be written as a general function of p and its gradients, constrained by symmetries and the 
necessity that the net energy decreases with time. Indeed, this was the approach Limkumnerd originally tooli^, but 
the complexity of the resulting theory and the multiplicity of terms allowed by symmetry led them to specialized^ to 
a particular choi ce mo tivated by the Peach-Kochlcr force — leading to the equation of motion previously developed 
by Acharya et alW^. 

The assumption that the net continuum dislocation density determines the evolution, however, is an uncontrollecP3 
and probably invalid assumption. fRef. 1671 have argued that the chaotic motion of dislocations may lead to a statistical 
ensemble that could allow a systematic theory of this type to be justified, but consensus has not been reached on 
whether this will indeed be possible.) The situation is less analogous to deriving the Navier-Stokes equation (where 
local equilibrium at the viscous length is sensible) than to deriving theories of eddy viscosity in fully developed 
turbulence (where unavoidable uncontrolled approximations are needed to subsume swirls on smaller scales into an 
effective viscosity of the coarse-grained system) . Important features of how dislocations are arranged in a local region 
will not be determined by the net Burgers vector density, and extra state variables embodying their effects are needed. 
In the context of dislocation dynamics, these state variables are usually added as SSDs and yield surfaces - although 
far more complex memory effects could in principle be envisioned. 

Let us write as the microscopic dislocation density (the sum of line-i5 functions along individual dislocations, 
as in Eq. ([I]) and following equations). For the microscopic density, allowing both glide and climb, the dislocation 
current is directly given by the velocity t;°(x) of the individual dislocation passing through x (see Eq. 36): 

(46) 



7° 



Let be the microscopic quantity coarse-grained density over a length-scale 



^;^(x)= J rf3yi^0(x + y)z/;^(y). 



where is a smoothing or blurring function. Typically, we use a normal or Gaussian distribution p^ 



(47) 



(48) 
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For our purposes, we can define the SSD density as the difference between the coarse-grained density and the micro- 
scopic density!^ 

p^^^(x) = /(x) - p^(x) = pO(x) - J dVp"(x + y)ii'^(y). (49) 

(Ref. i69l calls this quantity the dislocation fluctuation tensor field.) 

First, we address the question of SSDs, which we do not include in our simulations. In the paslP^l, we have argued 
that they do not contribute to the long-range stresses that drive the formation of the cell walls, and that the successful 
generation of cellular structures in our simplified model suggests that they are not crucial. Here we go further, and 
suggest that their density may be small on the relevant length-scales for cell-wall formation, and also that in a theory 
(like ours) with scale-invariant structures it would not be consistent to add them separately. 

What is the dependence of the SSD density on the coarse-graining scale? Clearly p'^ contains all dislocations; 
clearly for a bent single crystal of size L, contains only those dislocations necessary to mediate the rotation across 
the crystal (usually a tiny fraction of the total density of dislocations). As E increases past the distance between 
dislocations, cancelling pairs of Burgers vectors through the same grid face will leave the GNDs and join the SSDs. If 
the dislocation densities were smoothly varying, as is often envisioned on long length scales, the SSD density would 
be roughly independent of E except on microscopic scales. But, for a cellular structure with gross inhomogeneities in 
dislocation density, the SSD density on the mesoscale may be much lower than that on the macroscale . Ver y tangibly, 
if alternating cell walls separated by £ have opposite misorientations (as is quite commonly observecP^'^ , then the 
SSD density for E > £ will include most of the dislocations incorporated into these cell walls, while for S < £ the cell 
walls will be viewed as geometrically necessary. 

How does the GND density within the cell walls compare with the total dislocation density for a typical material? Is 
it possible that the GNDs dominate over SSDs in the regime where these cell wall patterns form? Recent simulations 
clearly suggest (see Ref. [63] [Figure 5]) that the distinction between GNDs and SSDs is not clear at the length scale 
of a micron, and with reasonable definitions GNDs dominate by at least an order of magnitude over the residual 
average SSD density. But what about the experiments? While more experiments are necessary to clarify this issue, 
the existing evidence supports that at mesoscales, SSDs at least are not necessarily dominant. In particular, Ref. 
observes that cell boundary structures exhibit Davdav/b = C where Dav is the average wall spacing and 9av is the 
average misorientation angle with C ~ 650 for 'geometrically necessary' boundaries (GNBs) and C '-^ 80 for 'incidental 
dislocation' boundaries (IDBs). The resulting dislocation density should scale as 

1 9av C 

where h is the average spacing between GNDs in the wall^SI. There are some estimates available from the literature. 
Referenc^^ tells us for pure aluminum that Dav is often observed to be Dav = 1 — 5/im which leads to roughly 
Pgnd ^ lO^-^ X (2.6 — 65)/m^ and Pq^jj one order of magnitude smaller. Similar estimates in Ref. [70] give Pq^d — 
lO^'* — 6 X 10^^/m^ for aluminum at von Mises strains of e = 0.2 and 0.6 respectively. The larger von Mises strains, the 
higher dislocation density. Typically, in highly deformed aluminum (e ~ 2.7), the total dislocation density is roughly 
10^^/m^ (see Ref. 150]) . While SSDs within a cell boundary may exist, it is clearly far from true that SSDs dominate 
the dynamics in these experiments. 

These TEM analyses of cell boundary sizes and misorientations have a misorientation cutoff 6*0 ~ 2'^^, they analyze 
the cell boundaries using a single typical length scale Dav Our model behavior is formally much closer to the fractal 
scaling analysis that Ref. [7] used. How does one identify a cutoff in a theory exhibiting scale invariance (i.e., with no 
natural length scale)? Clearly our simulations are cut off at the numerical grid spacing, and the scale invariant theory 
applies after a few grid spacings. Similarly, if the real materials are described by a scale-invariant morphology (still 
an open question), the cutoff to the scale invariant regime will be where the granularity of the dislocations becomes 
important - the dislocation spacing, or perhaps the annihilation length. This is precisely the length scale at which the 
dislocations are individually resolved - at which there is no separate populations of SSDs and GNDs. Thus ignoring 
SSDs in our theory is at least self-consistent. 

So, not only are the SSDs unimportant for the long-range stresses and appear unnecessary for our (presumably 
successful) modeling of the formation of cell walls, but they also may be rare on the sub-cellular coarse-graining scale 
we use in our modeling, and it makes sense in our mesoscale theory for us to omit their effects. 

The likelihood that we do not need to incorporate explicit SSDs in our equations of motion does not mean that our 



equations are correct. The microscopic equation of motion, Eq. (46 1 naively looks the same as our 'single- velocity' 



equation of motion we use {e.g., Eq. 36). But, as derived in Ref. i25i the coarse-graining procedure (Eq. |47[) leads to 
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a correction term to the single-velocity equations: 



= «fpf.,+ir,- (51) 

Acharya interpret^^Il this correction term as the strain rate due to SSD^^mH^ and later BeaudoirP^l and otherJ^ 
then use traditional crystal plasticity SSD evolution laws for it. Their GNDs thus move according to the same single- 
velocity laws as ours do, supplemented by SSDs that evolve by crystal plasticity (and thereby contribute changes to 
the GND density). This is entirely appropriate for scales large compared to the cellular structures, where most of the 
dislocations are indeed SSDs. 

Although we argue that SSDs are largely absent at the nanoscale where we are using our continuum theory, this does 
not mean the single-velocity form of our equations of motion can be trusted. Unlike fluid mixtures, where momentum 
conservation and Galilean invariance lead to a shared mean velocity after a few collision times, the microscopic 
dislocations are subject to different resolved shear stresses and are mobile along different glide planes, so neighboring 
dislocations may well move in a variety of directiond^. If so, the microscopic velocity will fluctuate in concert with 
the microscopic Burgers vector density on microscopic scales, and the correction will be large. Hence Acharya's 
correction term also incorporates multiple velocities for the GND density. Our single- velocity approximation (e.i?., 
Eq. 36 1 must be viewed as a physically allowed equation of motion, but a second uncontrolled approximation - the 
general evolution law for the coarse-grained system will be more complex. 

Let us be perfectly clear that our arguments, compelling on scales small compared to the mesoscale cellular struc- 
tures, should not be viewed as a critique of the use of SSDs on larger scales. Much of our understanding of yield stress 
and work hardening revolves around the macroscopic dislocation density, which perforce are due to SSDs (since they 
dominate on macroscopic scales). We also admire the work of Beaudoin, Acharya, and others which supplements the 



GND equations we both study with crystal plasticity rules for the SSDs motivated by Eq. ( 51 ). Surely on macroscales 
the SSDs dominate the deformation, and using a single-velocity law for the GNDs is better than ignoring them alto- 
gether, and we have no particular reason to believe that the contribution of multiple GND velocities in the evolution 
laws through will be significant or dominant. 

IV. RESULTS 
A. Two and three dimensional simulations 



We perform simulations in 2D and 3D for the dislocation dynamics of Eq. ( 15 1 and Eq. ( 18 1, with dynamical currents 



defined by CGD (Eq. |36|), GOD-MDP (Eq. |40|), and GOD-LVP (Eq. |44|. We numerically observe that simulations 
of Eqs. (15 1, (181 lead to the same results statistically (i.e., the numerica time step approximations leave the physics 
invariantJT^ We therefore focus our presentation on the results of Eq. (18 1, where the evolving field variable /3p is 
unconstrained. Our CGD and GOD-MDP models have been quite extensively simulated in one and two dimensions 
and relevant results can be found in Refs. ,22. .26^ and i73i In this paper, we concentrate on periodic grids of spatial 
extent L in both twcP3 and three dimensions. The numerical approach we use is a second-order central upwind 
scheme designed for Hamilton-Jacobi equation^^ using finite differences. This method is quite efficient in capturing 
(5— shock singular structures'^, even though it is flexible enough to allow for the use of approximate solvers near the 
singularities. 

Our numerical simulations show a close analogy to those of turbulent flows^. As in three-dime nsiona l turbulence, 
defect structures lead to intermittent transfer of morphology to short length scales. As conjecturecP^'^ fo r the Euler 
equations or the inviscid limit of Navier-Stokes equations, our simulations develop singularities in finite time^^'^^'. Here 
these singularities are 5-shocks representing grain-boundary-like structures emerging from the mutual interactions 
among mobile dislocation^^. In analogy with turbulence, where the viscosity serves to smooth out the vortex- 
stretching singularities of the Euler equations, we have explored the effects of adding an artificial viscosity term to 
our equations of motion'^ . In the presence of artificial viscosity, our simulations exhibit nice numerical convergence 
in all dimension^^. However, in the limit of vanishing viscosity, the solutions of our dynamics continue to depend 
on the lattice cutoff in higher dimensions, (our simulations only exhibit numerical convergence in one dimension). 
Actually, the fact that the physical system is cut off by the atomic scale leads to the conjecture that our equations 
are in some sense non-renormalizable in the ultraviolet. These issues are discussed in detail in Refs. l49l and [77l See 
also Ref. [75] for global existence and uniqueness results from an alternative regularization for this type of equations; 
it is not known whether these alternative regularizations will continue to exhibit the fractal scaling we observe. 
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FIG. 4: (Color online) Complex dislocation structures in two dimensions (1024^) for the relaxed states of an initially 
random distortion. Top: Dislocation climb is allowed; Middle: Glide only using a mobile dislocation population; Bottom: Glide 
only using a local vacancy pressure. Left: Net GND density \q\ plotted linearly in density with dark regions a factor ~ lO** 
more dense than the lightest visible regions, (a) When climb is allowed, the resulting morphologies are sharp, regular, and 
close to the system scale, (c) When climb is forbidden using a mobile dislocation population, there is a hierarchy of walls on 
a variety of length scales, getting weaker on finer length scales, (e) When climb is removed using a local vacancy pressure, the 
resulting morphologies are as sharp as those (a) allowing climb. Right: Corresponding local crystalline orientation maps, with 
the three components of the orientation vector A linearly mapped onto a vector of RGB values. Notice the fuzzier cell walls 
(c) and (d) suggests a larger fractal dimension. 
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FIG. 5: (Color online) Complex dislocation structures in three dimensions (128^) for the relaxed states of an initially 
random distortion. Notice these textured views on the surface of simulation cubes. Top: Dislocation climb is allowed; Middle: 
Glide only using a mobile dislocation population; Bottom: Glide only using a local vacancy pressure. Left: Net GND density | q\ 
plotted linearly in density with dark regions a factor ~ 10^ more dense than the lightest visible regions. The cellular structures 
in (a), (c), and (e) seem similarly fuzzy; our theory in three dimensions generates fractal cell walls. Right: Corresponding local 
crystalline maps, with the three components of the orientation vector A linearly mapped onto a vector of RGB values. 
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FIG. 6: (Color online) The elastic free energy decreases to zero as a power law in time in both two and three 
dimensions. In both (a) and (b), we show that the free energy T decays monotonically in time, and goes to zero as a power 
law for CGD, GOD-MDP, and GOD-LVP simulations, as the system relaxes in the absence of external strain. 



In the vanishing viscosity limit, our simulations exhibit fractal structure down to the smallest scales. When varying 
the system size continuously, the solutions of our dynamics exhibit a convergent set of correlation functions of the 
various order parameter fields, which are used to characterize the emergent self-similarity. This statistical convergence 
is numerically tested in |D 1[ 

In both two and three dimensional simulations, we relax the deformed system with and without dislocation climb 
in the absence of external loading. Here, the initial plastic distortion field /Jp is still a Gaussian random field with 
correlation length scale \/2L/b ^ 0.28L and initial amplitude /3o = 1. (In our earlier work^^, we described this 
length as L/5, using a non-standard definition of correlation length scale; see D2 ) These random initial conditions 
are explained in D2 In 2D, Figure |4] shows that CGD and GOD-LVP simulations (top and bottom) exhibit much 
sharper, flatter boundaries than GOD-MDP (middle). This difference is quantitatively described by the large shift 
in the static critical exponent rj in 2D for both CGD and GOD-LVP. In our earlier worlP^, we announced this 
difference as providing a sharp distinction between high-temperature, non- fractal grain boundaries (for CGD), and 
low-temperature, fractal cell wall structures (for GOD-MDP). This appealing message did not survive the transition 
to 3D; Figure [5] shows basically indistinguishable complex cellular structures, for all three types of dynamics. Indeed, 
Table H] shows only a small change in critical exponents, among CGD, GOD-MDP, and GOD-LVP. During both two 
and three dimensional relaxations, their appropriate free energies monotonically decay to zero as shown in Fig. [6] 



B. Self-similarity and initial conditions 

Self-similar structures, as emergent collective phenomena, have been studied in mesoscale crystal^^H^ human-scale 
social networl|2SI, and the astronomical-scale uniyers^^. In some model^Sfll, the self-similarity comes from scale-free 
initial conditions with a power-law spectrurrPi^. In our CDD model, our simulations start from a random plastic 
distortion with a Gaussian distribution characterized by a single length scale. The scale-free dislocation structure 
spontaneously emerges as a result of the deterministic dynamics. 

Our Gaussian random initial condition is analogous to hitting a bulk material randomly with a hammer. The ham- 
mer head (the dent size scale) corresponds to the correlated length. We need to generate inhomogeneous deformations 
like random dents, because our theory is deterministic and hence uniform initial conditions under uniform loading 
will not develop patterns. 

We have considered alternatives to our imposition of Gaussian random deformation fields as initial conditions, (a) As 
an alternative to random initial deformations, we could have imposed a more regular (albeit nonuniform) deformation 
- starting with our material bent into a sinusoidal arc, and then letting it relax. Such simulations produce more 
symmetric versions of the fractal patterns we see; indeed, our Gaussian random initial deformations have correlation 
lengths 'hammer size' comparable to the system size, so our starting deformations are almost sinusoidal (although 



different components have different phases) - see D 2 (b) To explore the effects of multiple uncorrelated random 
domains (multiple small dents), we reduce the Gaussian correlation length as shown in Fig. [jj We find that the 
initial-scale deformation determines the maximal cutoff for the fractal correlations in our model. In other systems 
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FIG. 7: (Color online) Relaxation with various initial length scales in two dimensions. GNDs are not allowed to climb 
due to the constraint of a mobile dislocation population in these simulations, (a), (b), and (c) are the net GND density map 
\g\, the net plastic distortion 1/3^ \ (the warmer color indicating the larger distortion), and the crystalline orientation map in a 
fully-relaxed state evolved from an initial random plastic distortion with correlated length scale 0.07L. They are compared to 
the same sequence of plots, (d), (e), and (f), which are in the relaxed state with the initial length scale 0.211/ three times as 
long. Notice the features with the lo ngest wave length reflecting the initial distortion length scales, (g), (h), and (i) are the 
scalar forms (discussed in Sec. IV C I of correlation functions of the GND density p, the intrinsic plastic distortion /?^'\ and 
the crystalline orientation A for well-relaxed states with initial length scales varying from 0.071/ to 0.28L. They exhibit power 
laws independent of the initial length scales, with cutoffs set by the initial lengths. (The scaling relation among their critical 
exponents will be discussed in Sec. [v|) 



(such as two-dimensional turbulence) one can observe an 'inverse cascade' with fractal structures propagating to long 
length scales; we observe no evidence of these here, (c) As an alternative to imposing an initial plastic deformation 
field and then relaxing, we have explored deforming the material slowly and continuously in time. Our pr elim inary 
'slow hammering' explorations turn the Gaussian initial conditions /Sp" into a source term, modifying Eq. ( 18 ) with 
an additional term to give dtPfj = Jij + Pf^ /t. Our early explorations suggest that slow hammering simulations will 
be qualitatively compatible with the relaxation of an initial rapid hammering. In this paper, to avoid the introduction 
of the hammering time scale r, we focus on the (admittedly less physically motivated) relaxation behavior. 

In real materials, initial grain boundaries, impurities, or sample sizes, can be viewed as analogies to our initial dents 
— explaining the observation of dislocation cellular structures both in single crystals and polycrystalline materials. 
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Figure [t] shows relaxation without dislocation climb (due to the constraint of a mobile dislocation population) 
at various initial length scales in 2D. From Fig. [7]ja) to (f), the net GND density, the net plastic distortion, and 
the crystalline orientation map, measured at two well-relaxed states evolved from different random distortions, all 
show fuzzy fractal structures, distinguished only by their longest-length-scale features that originate from the initial 
conditions. In Fig. [Tf^g), (h), and (i), the correlation functions of the GND density p, the intrinsic plastic distortion 
/3P'\ and the crystalline orientation A are applied to characterize the emergent self-similarity, as discussed in the 
following section [IV C[ They all exhibit the same power law, albeit with different cutoffs due to the initial conditions. 



C. Correlation functions 



Hierarchical dislocation structures have been observed both experimentalljM^ and in our simulation^^H Early 
work analyzed experimental cellular structures using the fractal box counting methocP or by separating the systems 
into cells and analyzing their sizes and misorientation^^lHlIl Jn our previous publication, we analyzed our simulated 
dislocation patterns using these two methods, and showed broad agreement with these experimental analyse^^Ii, 
fact, lack of the measurements of physical order parameters leads to incomplete characterization of the emergent 
self-similaritjff^. We will not pursue these methods here. 

In our view, the emergent self-similarity should best be exhibited by the correlation functions of the order parameter 
fields, such as the GND density p, the plastic distortion /3p, and the crystalline orientation vector A. Here we focus 
on scalar invariants of the various tensor correlation functions. 



For the vector correlation function C/^(x) (Eq. 52 ), only the sum C^(x) is a scalar invariant under three dimensional 
rotations. For the tensor fields p and ^p, their two-point correlation functions are measured in terms of a complete 
set of three independent scalar invariants, which are indicated by 'tot' ( total), ' per' (per mutation), and 'tr' (trace). 



In searching for the explanation of the lack of scaling22 for /?p (see Sec. IV C 3 and El), we checked whether these 



independent invariants might scale independently. In fact, most of them share a single underlying critical exponent, 
except for the trace-type scalar invariant of the correlation function of /3P'\ which go to a constant in well- relaxed 



states, as discussed in Sec. IV A2 
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FIG. 8: (Color online) Correlation functions of A in both two and three dimensions. In (a) and (b), red, blue, and 
green lines indicate CGD, GOD-MDP, and GOD-LVP simulations, respectively. Left: Correlation functions of A are measured 
in relaxed, unstrained 1024^ systems; Right: These correlation functions are measured in relaxed, unstrained 128^ systems. All 
dashed lines show estimated power laws quoted in Table [I] 



1. Correlation function of crystalline orientation field 



As dislocations self-organize themselves into complex structures, the relative differences of the crystalline orienta- 
tions are correlated over a long length scale. 

For a vector field, like the crystalline orientation A, the natural two-point correlation function is 

C,^(x) = ((A,(x)-A,(0))(A,(x)-A,(0))) 

= 2(A,Aj) -2(A,(x)A,(0)). (52) 
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FIG. 9: (Color online) Correlation functions of q in both two and three dimensions. Left: (a) is measured in relaxed, 
unstrained 1024^ systems; Right: (b) is measured in relaxed, unstrained 128^ systems. All dashed lines show estimated power 
laws quoted in Table |l] Notice all three scalar forms of the correlation functions of GND density share the same power law. 



Note that we correlate changes in A between two points. Just as for the height-height correlation function in surface 
growtfP^, adding a constant to A(x) (rotating the sample) leads to an equivalent configuration, so only differences in 
rotations can be meaningfully correlated. 
It can be also described in Fourier space 

C^(k) = 2(A,A,)(2^)35(k) - |,A,(k)A,(-k). (53) 

In an isotropic medium, we study the scalar invariant formed from Cfj 

C^(x) = C,^(x) = 2(A2) - 2(A,(x)A,(0)). (54) 

Figure [s] shows the correlation functions of crystalline orientations in both 1024^ and 128^ simulations. The large 
shift in critical exponents seen in 2D (Fig. [Sja)) for both CGD and GOD-LVP is not observed in the fully three 
dimensional simulations (Fig. |8][b)). 



2. Correlation function of GND density field 



As GNDs evolve into (5-shock singularities, the critical fluctuations of the GND density can be measured by the 
two-point correlation function Clx) of the GND density, which decays as the separating distance between two sites 
increases. The complete set of rotational invariants of the correlation function of p includes three scalar forms 

CL(x) = (p,,(x)p„-(0)), (55) 
Cp^,,(x) - (p,,(x)p,,(0)), (56) 
C(x) = (p,,(x)p„-(0)). (57) 

Figure [9] shows all the correlation functions of GND density in both 1024^ and 128"^ simulations. These three 
scalar forms of the correlation functions of p exhibit the same critical exponent ?/, as listed in Table |Tj Similar to the 
measurements of C"^, the large shift in critical exponents seen in 2D (Fig.|9];a)) for both CGD and GOD-LVP is not 
observed in the fully three dimensional simulations (Fig. [9j^b)). 



3. Correlation function of plastic distortion field 



The plastic distortion /3p is a mixture of both the divergence- free P^'^ and the curl-free (3^'^. Figure 10 shows that 
/3P does not appear to be scale invariant, as observed in our earlier work^^. It is crucial to study the correlations of 
the two physical fields, fi^'^ and /3P'^, separately. 
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Similarly to the crystalline orientation A, we correlate the differences between p'^-^ at neighboring points. The 
complete set of scalar invariants of correlation functions of (3^'^ thus includes the three scalar forms 



cf:/(x) 



3PJ 



P,I/ 



C"(x) = (x)-/3r(0))(/35^^(x)-/35^\0))) 



pj/ 



3P,I/ 



3P,I/ 



= 2i 



3P,IoP,I 



P,i^ 



/3P:0-2(/?r(x)/3]:;^(O)); 



(58) 



(59) 



(60) 



where an overall minus sign is added to C^^^ so as to yield a positive measure. 

the correlation functions of the intrinsic plastic distortion (3^'^ in both 1024^ and 128^ simulations exhibit 
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In Fig. 

a critical exponent cr'. Thes e measured critical exponents are shown in Table [l] We discuss the less physically relevant 



case of /3P^" in E 1 Fig. 17 



V. SCALING THEORY 

The emergent self-similar dislocation morphologies are characterized by the rotational invariants of correlation 
functions of physical observables, such as the GND density p, the crystalline orientation A, and the intrinsic plastic 
distortion /3P'^. Here we derive the relations expected between these correlation functions, and show that their critical 
exponents collapse into a single underlying one through a generic scaling theory. 

In our model, the initial elastic stresses are relaxed via dislocation motion, leading to the formation of cellular 
structures. In the limit of slow imposed deformations, the elastic stress goes to zero in our model. We will use the 
absence of external stress to simplify our correlation function relations. (Some relations can be valid regardless of 
the existence of residual stress.) Those relations that hold only in stress-free states will be labeled 'sf; they will be 
applicable in analyzing experiments only insofar as residual stresses are small. 




— CGD 

— GOD-MDP 

— GOD-LVP 



0.002 



R 



0.2 



FIG. 10: (Color online) Correlation functions of in two dimensions. Red, blue, and green lines indicate CGD, 
GOD-MDP, and GOD-LVP simulations, respectively. None of these curves shows a convincing power law. 
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— tot(CGD), o' = 1.1 

— per (CGD), o' = 1.15 

— tot (GOD-MDP), a' = 1.45 

— per (GOD-MDP), o' = 1.45 

— tot(GOD-LVP), o' = 1.1 
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0.01 R [L] 0.1 




R [L] 



FIG. 11: (Color online) Correlation functions of j3'^'^ in both two and three dimensions. In (a) and (b), the correlation 
functions of the intrinsic part of plastic distortion field are shown. Left: (a) is measured in relaxed, unstrained 1024^ systems; 
Right: (b) is measured in in relaxed, unstrained 128^ systems. All dashed lines show estimated power laws quoted in Table III 

Notice that we omit the correlation functions of Cj^ , which are independent of distance, and unrelated to the emergent 
self-similarity, as shown in Sec. |V A2[ 



A. Relations between correlation functions 

1. C andC^ 



For a stress- free state, we thus ignore the elastic strain term in Eq. ( 14 ) and write in Fourier space 



Pij{k) = -ikjAi{k) + iSijkkAk{k). 



(61) 



First, we can substitute Eq. (61 ) into the Fourier-transformed form of the correlation function Eq. (55) 



Ctot{k) = ^ { ~ikjAi(k) + i6^jkkAk{k) ] [ ikjAi{-k) - i%fc„A,„(-k) 



■u 1 

V 



{S,jk^ + hkj)A,{k)A,i-k). 



Multiplying both sides of Eq. (53) by {Sijk + kikj) gives 



{5,,e + hk,)C^A\^^ 



V 



iS,je + hkj)A,{k)A,{-k). 



Comparing Eq. (63) and Eq. (62), we may write C^^j in terms of C/^ as 



cu^)tl-hs,,k^ + k,k,)c^^ik). 



Second, we can substitute Eq. (61 1 into the Fourier-transformed form of the correlation function Eq. (56 1 



C;,,(k)l^|,fc,fc,A,(k)A,(-k). 



Multiplying both sides of Eq. (53) by kikj and comparing with Eq. (65) gives 



^perO^) — kikjC:(^{k). 



Finally, we substitute Eq. (61) into the Fourier-transformed form of the correlation function Eq. (57) 



C(k)=^ifc,fc,A,(k)A,(-k). 



(62) 



(63) 



(64) 



(65) 



(66) 



(67) 



22 



Repeating the same procedure of deriving C^g^, we write Cf^ in terms of C{j as 

Ct(k) H -2fc,fc/^(k). 



Through an inverse Fourier transform, we convert Eq. (64 1, Eq. ([66^, and Eq. (68 1 back to real space to find 



(68) 

(69) 

(70) 
(71) 



2. C^"' andC^ 



The intrinsic part of the plastic distortion field is directly related to the GND density field. In stress-free states, 
the crystalline orientation vector can fully describe the GND density. We thus can connect C^^' to C^. 
First, substituting (3fj^ = —isnmkiPmj/k'^ into the Fourier-transformed form of Eq. (58) gives 



C;(k) = 2(/3f//?P'^)(2^)^5(k) - U-^eum'^Pm,{\^)) ('^e..,^p,,(-k) 



3P,IoP,I\ 



2(/3S'^/3S'^)(27r)^<5(k)-^(^lp„,(k)p,„,(-k) 



(72) 



During this derivation, some terms vanish due to the geometrical constraint on p, Eq. ([6|. Multiplying — fc^/2 on both 
sides of Eq. (72) and applying the Fourier-transformed form of Eq. (55) gives 



cl7(k)=cL(k). 



(73) 



In stress-free states, we can substitute Eq. (64) into Eq. (73) 



which is rewritten after multiplying — 2/fc^ on both sides 



cr(k)i^'c^k) + ^c^ 



Second, substituting /3fj^ = —ieumkiPmj/k'^ into the Fourier-transformed form of Eq. (59) gives 

^per(k) 



-2{/3f.'pl;'){2nfS(k) + ^ U^e,^^^p,„,(k)j (ze,,,^p,,(-k) 



2{/3f:'l3^l'){27rrS{k} - ^fc,fc,p™,(k)p™(- 



+ ^2 ^tot (k) ^2 ^tr (k) , 



(74) 



(75) 



(76) 



where we skip straightforward but tedious expansions and the geometrical constraint on p, Eq. ([6|. Notice that this 
relation is correct even in the presence of stress. 



In stress-free states, we substitute Eqs. (61), (64), (68) into Eq. (76), and ignore the constant zero wavelength term 



ikjA^{k) + i(5,„jfcfcAfe(k) zfciA„(-k) - i(5„i/c„A„(-k) 



sf 



Vk^ 

-j^ik'^^ij + kikj)C^{k) + —kikjC^ik) 



CtAk). 



(77) 
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Finally, substituting /3f^^ = —ieumkiPmj/k'^ into the Fourier-transformed form of Eq. (60) gives 



2 ~ 2 ~ 



(78) 



valid in the presence of stress. Here we repeat a similar procedure as was used to derive in Eq. ( 76 ) 
In stress-free states, we substitute Eqs. (61 1, (64), (66) into Eq. (78) 



Ct%) 2{/3f^l3^f}(2n)'S(k) + ^{kH,, + hk,)C^^{\.) - ^hk,C^{k) 



1 

fc2' 



2 ki k j 
' Vk^ 



«fciA,„(k) -I- i(5„iifcfcAfe(k) ifcjA,„(-k) - i(5„y/c„A„(-k) 



2{l3f^'m')i27rf6{k), 



which is a trivial constant in space. 



Through an inverse Fourier transform, Eqs. (75), (77), and (79) can be converted back to real space, giving 



R3 3 )C,^(x'), 



CrV)=2/d3x'/3P'i(x')/3P;^(x') = 2(;9f/?P;i), 



where R = x' — x. According to Eqs. (75) and (77), we can extract a relation 



(x) - 2Ct^, (x) + 2C^(x) const 



(79) 

(80) 
(81) 
(82) 

(83) 



TABLE I: Critical exponents for correlation functions at stress-free states. (C.F. and S.T. represent 'Correlation 
Functions' and 'Scaling Theory', respectively.) 



C.F. 



S.T. 



Climb&Glide 



Simulations 
Glide Only (MDP) 



LVP Glide Only (LVP) 











2D(1024^) 


3D(128^) 


2D(1024^) 


3D(128^) 


20(1024^) 


30(128=*) 


'-'tot 




V 




0.80 ±0.30 


0.55 ±0.05 


0.45 ±0.25 


0.60 ±0.20 


0.80 ±0.30 


0.55 ±0.05 


c 




V 




0.80 ±0.20 


0.55 ±0.05 


0.45 ±0.20 


0.60 ±0.20 


0.70 ±0.30 


0.50 ±0.05 


C'tr 




V 




0.80 ±0.20 


0.55 ±0.05 


0.45 ±0.20 


0.60 ±0.10 


0.70 ±0.30 


0.45 ±0.05 




2 




V 


1.10 ±0.65 


1.45 ±0.25 


1.50 ±0.30 


1.35 ±0.25 


1.10 ±0.65 


1.50 ±0.25 


'-'tot 


2 




V 


1.10 ±0.60 


1.45 ±0.15 


1.45 ±0.25 


1.30 ±0.20 


1.10 ±0.60 


1.50 ±0.20 


^/3P.I 


2 




V 


1.15 ±0.45 


1.50 ±0.25 


1.45 ±0.25 


1.50 ±0.50 


1.20 ±0.45 


1.55 ±0.25 



We can convert Eq. ( 73 ) through an inverse Fourier transform 



C,(x) = i52C(x), 



or 



cL (x) 



27r ./ i? ' 



(84) 



(85) 



valid in the presence of residual stress. 
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B. Critical exponent relations 



When the self-similar dislocation structures emerge, the correlation functions of all physical quantities are expected 
to exhibit scale-free power laws. We consider the simplest possible scenario, where single variable scaling is present 
to reveal the minimal number of underlying critical exponents. 

First, we define the critical exponent rj as the power law describing the asymptotic decay of Cfot(x) ^ Ixj^**, one of 
the correlation functions for the GND density tensor (summed over components). If we rescale the spatial variable x 
by a factor b, the correlation function C is rcscaled by the power law as 



Similarly, the correlation function of the crystalline orientation field A is described by a power law, C'^(x) 
where a is its critical exponent. We repeat the rescaling by the same factor b 

C^(&x) = 6"C^(x). 



(86) 



(87) 



Since Cf^j can be written in terms of C^, Eq. (69 1, we rescale this relation by the same factor b 



cL(&x) 



sj 1 

2 



C^(6x) + - 



9, 



Substituting Eq. (87) into Eq. (88 1 gives 



/ icr-2 



^92C^(x) + ia,9,C,^(x) 



cL(x). 



Comparing with Eq. ( 86 1 gives a relation between a and rj 



(89) 
(90) 



We can repeat the same renormalization group procedure to analyze the critical exponents of the other two scalar 
forms of the correlation functions of the GND density field. Clearly, C^^^ and Cf,. share the same critical exponent f] 
with Cfot. 

Also, we can define the critical exponent a' as the power law describing the asymptotic growth of Cf^^ (x) ^ |x|°' , 
one of the correlation functions for the intrinsic part of the plastic distortion field. We can rescale the correlation 
function C^^' 



Cf„7(6x)=fe'^'C;(x) 



(91) 



We rescale the relation Eq. (84 1 by the same factor 6, and substitute Eq. (|91|) into it 

1 



CL(6x) = 



CL (b^) = b' 



CM- 



Comparing with Eq. ( 86 1 also gives a relation between cr' and rj 



a' = 2-7]. 

op, I . 

Since both Ct^j and C share the same critical exponent 2 — 



V: 



it is clear that C^^^ 



(92) 



(93) 



the other scalar form of the 



correlation functions of the intrinsic plastic distortion field, also shares this critical exponent, according to Eq. (83). 



Thus the correlation functions of three physical quantities (the GND density p, the crystalline orientation A, 
and the intrinsic plastic distortion f3^'^) all share the same underlying universal critical exponent i] for self-similar 
morphologies, in the case of zero residual stress, and still hold in the limit of slow imposed deformation. Table |T] 
verifies the existence of single underlying critical exponent in both two and three dimensional simulations for each 
type of dynamics. Imposed strain, studied in Ref. 1221 could in principle change rj, but the scaling relations derived 
here should still apply. The strain, of course, breaks the isotropic symmetry, allowing even more allowed correlation 
functions to be measured. 
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C. Coarse graining, correlation functions, and cutoffs 



Our dislocation density p, as discussed in Sec. |III[ is a coarse-grained average over some distance E - taking the 
discrete microscopic dislocations and yielding a continuum field expressing their flux in different directions. Our power 
laws and scaling will be cut off in some way at this coarse-graining scale. For our simulations, the correlation functions 
extend down to a few times the numerical grid spacing (depending on the numerical diffusion in the algorithm we use) . 
For experiments, the correlation functions will be cut off in ways that are determined by the instrumental resolution. 
Since the process of coarse-graining is at the heart of the renormalization-group methods we rely upon to explain the 
emergent scale invariance in our model, we make an initial exploration here of how coarse-graining by the Gaussian 
blur of Eq. ( 47 ) and Eq. ( 48 ) affects the — correlation function. 



Following Eq. (47), 



cL(x) = {pf,{^)pm) 



1 



1 



V (27rl]2) 



/ d'zp^,(y + z)e 



dVp?.(y + x + z')e-^"/(2E^). 



(94) 



By changing variables s = y -f- z and A = z' — z, we integrate out the variable z of Eq. (94) 
cCi^) = ^Z^^ J J d^sp%is)pl{s + A + ^)e- 



87r3/2 
1 



87r3/2E3 



-AV(4S^) 



(95) 



In our simulating system, the correlation functions of GND density can be described by a power-law (Eq. 86 1, 
Cfof(x) = gjxl^'', where g is a constant. Thus, Eq. (95) is 



Cfo:(x) = ^/rf^A|x-fA|-e--V(4.^). 



(96) 



This correlation function of the coarse-grained GND density at the given scale S is a power-law smeared by a Gaussian 
distribution. 

Since the scalar field of the coarse-grained correlation function is rotational invariant, we assume that x is aligned 
along the x axis, x = {x, 0, 0). Then we could evaluate the integral of Eq. (96 1 in cylindrical coordinates A = {X, r, 9) 



9 



87r3/2E3 
9 



27Trdr / dX|(x + X)2+r2r"/2e-(^'+'-')/(4S^) 



^^^S-"-! / dXe(^ +2.x)/(4S^)p(^ _ ^/2, {x + X) V(4S2)) . 



(97) 



We can rewrite this coarse-grained correlation functions Eq. ( 97 ) as a power-law multiplied by a scaling function 

cC{x,J:) ^ g\x\-^^{x/^), 



where the scaling function ^'(•) (Fig. 12) equals 

^ J —OO 



(98) 



(99) 



VI. CONCLUSION 

In our earlier workj^ 26 | 49 |^ have proposed a flexible framework of ODD to study complex mesoscale phenomena 
of collective dislocation motion. Traditionally, deterministic CDDs have missed the experimentally ubiquitous feature 
of cellular pattern formation. Our CDD models have made progress in that respect. In the beginning, we focused our 
efforts on describing coarse-grained dislocations that naturally develop dislocation cellular structures in ways that are 
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FIG. 12; Scaling function of the correlation function of coarse-grained GND density p^. We calculate the correlation 
function of the coarse-grained GND density at the given scale E. Theoretically, its scaling function remains a power-law at 
the small coarse-graining length scale, and flattens out to be 1 as the correlation length of the system is far larger than the 
coarse-graining scale. 

consistent with experimental observations of scale invariance and fractality, a target achieved in Ref. ,22, However, 
that paper studied only 2D, instead of the more realistic 3D. 

In this manuscript, we go further in many aspects of the theory extending the results of our previous work: 

We provide a derivation of our theory that explains the differences with traditional theories of plasticity. In addition 
to our previously studied climb-glide (CGD) and glide-only (GOD-MDP) models, we extend our construction in order 
to incorporate vacancies, and re-deriv^^ a different glide-only dynamics (GOD-LVP) which we show exhibits very 
similar behavior in 2D to our CGD model. It is worth mentioning that in this way, the GOD-LVP and the CGD 
dynamics become statistically similar in 2D, while the previously studied, less physical, GOD-MDP model provides 
rather different behavior in 2Tj^. 

We present 3D simulation results here for the first time, showing qualitatively different behavior from that of 2D. 
In 3D, all three types of dynamics - CGD, GOD-MDP and GOD-LVP - show similar non-trivial fractal patterns and 
scaling dimensions. Thus our 3D analysis shows that the flatter 'grain boundaries' we observe in the 2D simulations 
are not intrinsic to our dynamics, but are an artifact of the artificial z-independent initial conditions. Experimentally, 
grain boundaries are indeed flatter and cleaner than cell walls, and our theory no longer provides a new explanation 
for this distinction. We expect that the dislocation core energies left out of our model would flatten the walls, and that 
adding disorder or entanglement would prevent the low-temperature glide-only dynamics from flattening as much. 

We also fully describe, in a statistical sense, multiple correlation functions - the local orientation, the plastic 
distortion, the GND density - their symmetries and their mutual scaling relations. Correlation functions of important 
physical quantities are categorized and analytically shown to share one stress-free exponent. The anomaly in the 
correlation functions of which was left as a question in our previous publication'^, has been discussed and 
explained. All of these correlation functions and properties are verified with the numerical results of the dynamics 
that we extensively discussed. 

As discussed in Sec. |lj our model is an immensely simplified caricature of the deformation of real materials. How 
does it connect to reality? 

First, we show that a model for which the dynamics is driven only by elastic strain produces realistic cell wall 
structures even while ignoring slip systems, crystalline anisotropj^, pinning, junction formation, and SSDs. The fact 
that low-energy dislocation structures (LEDS) provides natural explanations for many properties of these structures 
has long been emphasized by Ref. |83l Intermittent flow, forest interactions, and pinning will in general impede access 
to low energy states. These real-world features, our model suggests, can be important for the morphology of the cell 
wall structures but are not the root cause of their formation nor of their evolution under stress (discussed in previous 
worlPl). 

One must note, however, that strain energy minimization does not provide the explanation for wall structures in 
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our model material. Indeed, there is an immense space of dislocation densities which make the strain energy zercPSl^ 
including many continuous densities. Our dynamics relaxes into a small subset of these allowed structures - it is the 
dynamics that leads to cell structure formation here, not purely the energy. In discrete dislocation simulations and 
real materials, the quantization of the Burgers vector leads to a weak logarithmic energetic preference for sharp walls. 
This — /ib/(47r(l — t^))6'log6' energy of low-angle grain boundaries yields a log 2 preference for one wall of angle 9 
rather than two walls of angle 9/2. This leads to a 'zipping' together of low angle grain boundaries. Since b — > in 
a continuum theory, this preference is missing from our model. Yet, we still find cell wall formation suggesting that 
such mechanisms are not central to cell wall formation. 

Second, how should we connect our fractal cell wall structures with those (fractal or non-fractal) seen in experiments? 
Many qualitatively different kinds of cellular structures are seen in experiments - variously termed cell block structures, 
mosaic structures, ordinary cellular structures, .... Ref. (35] recently categorized these structures into three types, 
and argue that the orientation of the stress with respect to the crystalline axes largely determines which morphology 
is exhibited. The cellular structures in our model, which ignores crystalline anisotropy, likely are the theoretical 
progenitors of all of these morphologies. In particular, Hansen's type 1 and type 3 structures incorporate both 
'geometrically necessary' and 'incidental dislocation' boundaries (GNBs and IDBs), while type 2 structures incorporate 
only the latter. Our simulations cannot distinguish between these two types, and indeed qualitatively look similar 
to Hansen's type 2 structures. One should note that the names of these boundaries are misleading - the 'incidental' 
boundaries do mediate geometrical rotations, with the type 2 boundaries at a given strain having similar average 
misorientations to the geometrically necessary boundaries of type 1 structureJ^SI (Figure 8). It is commonly asserted 
that the IDBs are formed by statistical trapping of stored dislocations; our model suggests that stochasticity is not 
necessary for their formation. 

Third, how is our model compatible with traditional plasticity, which focuses on the total density of dislocation lines? 
Our model evolves the net dislocation density, ignoring the geometrically unnecessary or statistically stored dislocations 
with cancelling Burgers vectors. These latter dislocations are important for yield stress and work hardening on 
macroscales, but are invisible to our theory (since they do not generate stress). Insofar as the cancellation of Burgers 
vectors on the macroscale is due to cell walls of opposing misorientations on the mesoscale, there needs to be no conflict 
here. Also our model remains agnostic about whether cell boundaries include significant components of geometrically 
unnecessary dislocations. However, our model does assume that the driving force for cell boundary formation is the 
motion of GNDs, as opposed to (for example) inhomogeneous flows of SSDs. 

There still remain many fascinating mesoscale experiments, such as dislocation avalanche^smH^ size-dependent 
hardness (smaller is stronger Jpsl^ and complex anisotropic loadin^^^lSSi^ ^j-^^t we hope to emulate. We intend in the 
future to include several relevant additional ingredients to our dynamics, such as vacancies (CI), impurities (C2|, 
immobile dislocations/SSDs and slip systems, to reflect real materials. 
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Appendix A: Physical quantities in terms of the plastic distortion tensor jS^ 

In an isotropic infinitely l arge medium, the local deformation u, the elastic distortion (5° and the internal long-range 
stress ct'"* can be expressecP^'^ in terms of the plastic distortion field /3p in Fourier space: 



S,(k) - iV,,,(k)/3P(k), 



Niki{k) = --^{kkSii + kiSik) - i 
^°.(k) - T,,fe;(k)^P (k), 



(l-!.)fc2 (l-iy)/c4' 



(1 - v)k^ 

ijmn 



ajf (k) = M,,™„(k)/3P„(k), 



A:2 



kjkji A A ^ 



fc2 



^:'hL k'l^ k J k^fy^ k^ 



(Al) 



(A2) 



(A3) 



All these expressions are valid for systems with periodic boundary conditions. 

According to the definition Eq. ( 12 ) of the crystalline orientation A, we can replace uj'^ with /3'^ and e" by using the 
elastic distortion tensor decomposition Eq. ( 10 ) 



(A4) 



Here the permutation factor acting on the symmetric elastic strain tensor gives zero. Hence we can express the 
crystalline orientation vector A in terms of (5^ by using Eq. (A2) 



Ai(k) = ■:j^eijk^-j^{kjks5kt + kjkt5ks-k'^5jsSkt) 



1 



(1 - l/)fc4 



2/fc2 



[sijtkjkg + Sijskjkt — fc^e^sf )/3^((k). 



(A5) 



Appendix B: Energy dissipation rate 
1. Free energy in Fourier space 

In the absence of external stress, the free energy J- is the elastic energy caused by the internal long-range stress 

T = Jd^^ i^lf = l^^x la.mneyi^r., (Bl) 

where the stress is (T^ — Cijmn^^tii with. Cijmn 

the stiffness tensor. 

Using the symmetry of Cijmn and ignoring large rotations, — {/3fj + /3|j)/2, we can rewrite the elastic energy 
in terms of (3° 
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Performing a Fourier transform on both (3f- and /3Pj„ simultaneously gives 



(27r)3 7 (27r)3 



:Cij„i„/3°j(k)/3^„(k ) 



(B3) 



Integrating out the spatial variable x leaves a (5— function (5(k + k') in Eq. (B3|. We hence integrate out the k-space 
variable k' 



d^k 1 ~ ~ 

7, CijmnPi^ (}^)Pmn {^'^) ■ 



Substituting Eq. ( A2 ) into Eq. ( B4 ) gives 

d^k 1 



(27r)3 2 
d^k 1 
(27r)3 2 



(27r)3 2 ^^™"^-''^ 

(k)T™„,t(-k))/3P,(k)/?P(-k) 
M,,,*(k)^P,(k)^P(-k), 



where we skip straightforward but tedious simplifications. 



When turning on the external stress, we repeat the same procedure used in Eq. (B3|, yielding 

d3k 



a,f(k)/3f,(-k). 



(27r)3 *J ^"''^^J^ 



(B4) 



(B5) 



(B6) 



2. Calculation of energy functional derivative with respect to the GND density q 



According to Eq. (171, the infinitesimal change of the variable 5q is given in terms of 5j3^ 

5Qijk = -g.jisdiiSPl^). 



(B7) 



rewritten in terms of 



/3P 



Substituting Eq. (B7) into Eq. (27) and applying integration by parts, the infinitesimal change of T is hence 



According to Eq. ( 24 ) , it suggests 



Comparing Ec^. ( B8 1 and Eq. ( B9 ) implies 



up to a total derivative which we ignore due to the use of periodic boundary conditions. 



(B8) 



(B9) 



(BIO) 



3. Derivation of energy dissipation rate 

We can apply variational methods to calculate the dissipation rate of the free energy. As is well known, the general 
elastic energy £ in a crystal can be expressed as £ = ^ J d^x (Jijelj, with e^^ the elastic strain. An infinitesimal change 
of £ is: 



(fx (JijSeij, 



(Bll) 
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where we use aijSe^j = C^jkieliSe^.j = Saijclj. 

So the infinitesimal change of the free energy Eq. ( 20 ) is 

We apply the relation e° = e — e^, where eP is the plastic strain and e is the total strain: 



(B12) 



(B13) 



Using the symmetry of aij and ignoring large rotations, — -^{diUj + djUi), we can rewrite the first term of 
Eq. (B13) as /d^x afj^d{dtUj). Integrating by parts yields /d^x (9i((5ujcr"*) — Sujdialf). We can convert the first 
volume mtegral to a surface integral, which vanishes for an infinitely large system. Hence 



J' 



ij 



(B14) 



The first term of Eq. (B14| is zero assuming instantaneous elastic relaxation due to the local force equilibrium 
condition, 



5F - 



d^x « + a-*)<5;9P. 
The free energy dissipation rate is thus 5F/5t for Sjif^ = ^^St, hence 



using the symmetry of a.ij and efj = ^{Pfj + /3jj) 



(B15) 



int , ^cxtN f^tj 



' dt 

d^x «+a-*)J,,. 



(B16) 



When dislocations are allowed to climb, substituting the CGD current Eq. (36) into Eq. (B16I implies that the free 
energy dissipation rate is strictly negative 



~dt 



D - 



(B17) 



When removing dislocation climb by considering the mobile dislocation population, we substitute Eq. (40) into 
Eq. (B16) to guarantee that the rate of the change of the free energy density is also the negative of a perfect square 



~dt 



(B18) 



Appendix C: Model Extensions: Adding vacancies and disorder to CDD 



1. Coupling vacancy diffusion to CDD 



In plastically deformed crystals at low temperature, dislocations usually move only in the glide plane because 
vacancy diffusion is almost frozen out. When temperature increases, vacancy diffusion leads to dislocation climb 
out of the glide plane. At intermediate temperatures, slow vacancy diffusion can enable local creep. The resulting 
dynamics should couple the vacancy and dislocation fields in non-trivial ways. Here we couple the vacancy diffusion 
to the dislocation motion in our CDD model. 
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We introduce an order parameter field c(x), indicating the vacancy concentration density at the point x. The free 
energy J" is thus expressed 

T = T^^^ + ^^''^ = 1 d^xQa,,4 + \a{c - c^f^ , (CI) 

where a is a positive material parameter related to the vacancy creation energy, and Cq is the overall equilibrium 
vacancy concentration density. 

Assuming that GNDs share the velocity v in an infinitesimal volume, we write the current J for GNDs 

The current trace Ju describes the rate of volume change, which acts as a source and sink of vacancies. The coupling 
dynamics for vacancies is thus given as 



dtc^jW^c + Ju, (C3) 



where 7 is a positive vacancy diffusion constant. 



The infinitesimal change of the free energy (Eq. CI ) is 



We apply Eq. ( |B15| ) and bT^^-'-jbc = a{c - cq) 

= y d^^ (^-<y^JSpf^ + a{c - co)6c 

The free energy dissipation rate is thus 5T /5t for Sjif^ — ^^St and 6c = ff (5c, hence 

/■ / dPl . . dc 



(C5) 



Substituting the cmrent J (Eq. \G2\ and Eq. ( C3 1 into Eq. ( C6 ) gives 

(f-x.{aij{vuQuij) - a(c- co)(7V^c + 
= - / (f'y.{{a.ij - a{c- CQ)6ij)Qmj)vu- d^xaj(yc)^, (C7) 



Ik 



where we integrate by parts by assuming an infinitely large system. 

If we choose the velocity w„ — ~ a{c — CQ)Sij') Qy^ij , [D is a positive material dependent constant and 1/\q\ is 

added for the same reasons as discussed in Sec. |II CT |, the free energy is guaranteed to decrease monotonically. The 



coupling dynamics for both GNDs and vacancies is thus 

I 5t/3fj = ifl (c^mn - a{c - Co)5rnn) QuninQmj , ^^^^ 

This dynamics gives us a clear picture of the underlying physical mechanism: the vacancies contribute an extra 
hydrostatic pressure p — —a{c — cq). 

2. Coupling disorder to CDD 

In real crystals, the presence of precipitates or impurities results in a force pinning nearby dislocations. We can 
mimic this effect by incorporating a spatially varying random potential field V^(x). 

In our CDD model, we can add the interaction energy between GNDs and random disorder into the free energy 
(Eq.|20l) 

F = Fe+Ti= j rf^^xQajf ~ + V{x)\b\\ , (C9) 
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where Fe indicates the elastic free energy corresponding to the integral of the first two terms, and !Fi indicates the 
interaction energy, the integral of the last term. 

An infinitesimal change of the free energy is written 



In an infinitely large system, Eq. (B15) gives 



5Te 



and Eq. (B8| implies 



6Tj 



(CIO) 



(Cll) 



d ^9iji 



(C12) 



Substituting Eq. ( |Cll| and Eq. (|C12|) into Eq. (|C10|) gives 

5 J- = 



(C13) 



where the effective stress field 



off — int I ^cxt 



5mnH<9;(y(x) 



lei 



By replacing <Tij with cr? in the equation of motion of either allowing climb (Eq. 36 1 or removing climb (Eqs 



and 44 ) , we achieve the new CDD model that models GNDs interacting with disorder 
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FIG. 13: (Color online) Statistical convergence of correlation functions of A, p and fi^'^ by varying lattice sizes 
in two dimensions. We compare correlation functions of relaxed glide-only states (GOD-MDP) at resolutions from 128^ to 
1024^ systems. Top: We see that the correlation functions in all cases exhibit similar power laws in (a), (b), and (c); Bottom: 
(d), (e), and (f) show a single underlying critical exponent which appears to converge with increasing resolution, where a is the 
grid spacing. The black dashed lines are guides to the eye. 
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FIG. 14: (Color online) Statistical convergence of correlation functions of A, p and fP'^ by varying the initial 
length scales in 1024^ simulations. We measure correlation functions of relaxed glide-only states (GOD-MDP) at initial 
correlated lengths from 0.07L to 0.28L. In (a), (b), and (c), the radial-length variable R is rescaled by their initial correlation 
lengths, and the corresponding correlation functions are divided by the same lengths to the exhibiting powers. They roughly 
collapse into the scaling laws. Notice that the power laws measured in the state with the initial correlated length 0.071/ get 
distorted due to the small outer cutoff. 



Appendix D: Details of the Simulations 
1. Finite size effects 

Although we suspect that our simulations don't have weak solution^^, we can show that these solutions converge 
statistically. We use two ways to exhibit the statistical convergence. 

When we continue to decrease the grid spacing to zero (the continuum limit), we show the statistical convergence of 



13 



correlation functions of p, A, and /J^'^, with a slow expected drift of apparent exponents with system size, see Fig 

We can also decrease the initial correlated length scales in a large two dimensional simulation. Since the emergent 
self-similar structures are always developed below the initial correlated lengths, as discussed in Sec. |IVB[ this is 
similar to decreasing the system size by reducing the initial correlated lengths. In Fig. [T4j the correlation functions 
of p, A, and /3P'^ collapse into a single scaling curve, using finite size scaling. 



2. Gaussian random initial conditions 



Gaussian random fields are extensively used in physical modelings to mimic stochastic fluctuations with a correlated 
length scale. In our simulations, we construct an initially random plastic distortion, a nine-component tensor field, 
where every component is an independent Gaussian random field sharing a underlying length scale. 

We define a Gaussian random field / with correlation length ctq by convolving white noise (^(x)^(x')) = 5(x — x') 
with a Gaussian of width a^. 

/(x) = j d3x'e(x')e-(^-^')'/'^o. (Dl) 

In Fourier space, this can be done as a multiplication: 

/(k) - e--o'='/4e(k). (D2) 

The square /(k)/(-k) = e-'^o^Va jj^plies that the correlation function (/(x)/(x')) = (27rcr§)-3/2e-(^-''')'/(2<^o). 
In our simulations, the initial plastic distortion tensor field fi^ is constructed in Fourier space 

^P(k) = e--"'=^/^4(k), (D3) 

where the white noise signal Q is characterized as (C(ij)(x)C(i j)(x')) = j)(5(x — x'), and in Fourier space 
(k)C(j (— k) = -^(i.j)- (We use (i,j) to indicate a component of the tensor field, to avoid the Einstein summa- 
tion rule.) The correlation function of each component of /3P'^ is thus expressed in Fourier space 

C = 2(/3[:;^)/3[:^))(2.)3<5(k)~|^,^,(k)^^^,(-k) 

= 2(/3P;^)/3P^))(27^)3<5(k)-2A(,,,)e--o^^/^ (D4) 
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FIG. 15: (Color online) Gaussian random initial conditions with the correlated length scale 0.281/ in two dimen- 
sions, (a) shows the initial net GND density map; (b) exhibits the correlation functions of p under various initial conditions, 
where we compare the Gaussian random field to both a sinusoidal wave and a single periodic superposition of Gaussian peaks. 
The kink arises due to the edges and corners of the square unit cell. 



where the Gaussian kernel width uq, as a, standard length scale, defines the correlation length of our simulation. (In 
our earlier work, we use a non-standard definition for the correlation length, so our co equals the old length scale 
times \/2.) 

According to Eq. ^ and Eq. (D3), we can express the initial GND density field p in Fourier space 



p,,(k) = -leame-"-'' /^fc,C™,(k). (D5) 
The scalar invariant C^g^ of the correlation function of p is thus expressed in Fourier space 

= ^Pij{\<L)pij{~'k) 

^ ^e-'^«'='/2(A:2<S„„-fc™fc„)C™,(k)C„,(-k). (D6) 

The resulting initial GND density is not Gaussian correlated, unlike the initial plastic distortion. Figure [T5] exhibits 
the initial GND density map due to the Gaussian random plastic distortions with the correlation length 0.28L, and 
its correlation function. We compare the latter to the correlation functions of both a sinusoidal wave and a single 
periodic superposition of Gaussian peaks. The similarity of the three curves shows that our Gaussian random initial 
condition at ctq ^ 0.28i approaches the largest effective correlation length possible for periodic boundary conditions. 



TABLE II: Critical exponents for correlation functions of strain-history-dependent fields at stress-free 
states. (C.F. and Exp. represent 'Correlation Functions' and 'Exponents', respectively.) 



C.F. 


Exp. 


Chmb&Glide 
20(1024^) 3D(128^) 


Simulations 
Glide Only (MDP) 
2D(1024^) 3D(128^) 


Glide Only (LVP) 
2D(1024^) 3D(128^) 


CP"''' 


r 
/ 

T 

t 

T 
n 

T 


0.65 ± 1.00 
0.70 ±0.95 
0.70 ±0.95 
1.90 ±0.10 


1.05 ±0.65 
1.10 ±0.60 
1.10 ±0.60 
1.85 ±0.15 


1.25 ±0.60 
1.95 ±0.05 
1.95 ±0.05 
1.95 ±0.05 


1.20 ±0.50 
1.75 ±0.15 
1.75 ±0.15 
1.90 ±0.10 


0.55 ± 1.10 
0.50 ± 1.15 
0.50 ± 1.15 
1.95 ±0.05 


1.05 ±0.65 
1.05 ±0.70 
1.05 ±0.70 
1.90 ±0.10 
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(a) 



(c) 



(e) 




(b) 



(d) 



(f) 




FIG. 16: (Color online) Strain-history-dependent fields /S^'^ and i/> in two dimensions for the relaxed states. Top: 
Dislocation climb is allowed; Middle: Glide-only using a mobile dislocation population; Bottom: Glide-only using a local 
vacancy pressure. Left: The strain-history-dependent plastic distortion (a), (c), and (e) exhibit patterns reminiscent 

of self-similar dislocation structures. Right: The strain-history-dependent plastic deformation lip], (b), (d), and (f) exhibit 
smooth patterns with a little distortion, which are not fractal. 
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(a) 




le-4^ 



M tot (CGD), X = 0.65 
GK)per(CGD), x' = 0.7 
M tr (CGD), X' = 0.7 
o-o tot (GOD-MDP), X = 1.25 
G-Oper (GOD-MDP), x' = 1.95 
M tr (GOD-MDP), x' = 1.95 
tot (GOD-LVP), X = 0.55 
oe per (GOD-LVP), x' = 0.5 
M tr (GOD-LVP), X' = 0.5 




0.002 



[L] 



0.2 



0.02 



/? [L] 0.2 



FIG. 17: (Color online) Correlation functions of fF'^ in both two and three dimensions. In both (a) and (b), the 
correlation functions of the strain-history-dependent part of the plastic distortion fi^''^ are shown. Left: (a) is measured in 



relaxed, unstrained 1024 systems; Right: 
estimated power laws quoted in Table [III 



(b) is measured in in relaxed, unstrained 128 systems. All dashed lines show 




0.002 



R [L] 




R [L] 



FIG. 18: (Color online) Correlation functions of xf> in both two and three dimensions. In (a) and (b), the correlation 
functions of the strain-history-dependent deformation tp are shown. Red, blue, green lines indicate CGD, GOD-MDP, and GOD- 
LVP, respectively. Left: (a) is measured in relaxed, unstrained 1024'^ systems; Right: (b) is measured in relaxed, unstrained 
128'^ systems. All dashed lines show estimated power laws quoted in Table [ll] 



Appendix E: Other correlation functions unrelated to static scaling theory 
1. Correlation functions of the strain-history-dependent plastic deformation and distortion fields 



The curl-free strain-history-dependent part of the plastic distortion field, as shown in Fig. 16 'a), (c), and (e). 



exhibits structures reminiscent of self-similar morphology. We correlate their differences at neighboring points 

(0))), 



(0))). 



'"tot 


(x) 








per 


(x) 








'"tr 


(x) 


= ((/3&"(x) 







(El) 
(E2) 
(E3) 



Consider also the deformation field ij} (shown in Fig. [l6|b), (d), and (f)) of Eq. (19) whose gradient gives the 
strain- history-dependent plastic deformation /3P'^. Similarly to the crystalline orientation A, we correlate differences 
of i/j. The unique rotational invariant of its two-point correlation functions is written 

C^(x) = 2(V;2)-2(7/;,(x)V',(G)). (E4) 
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In Fig. 17 the correlation functions of the strain-history-dependent plastic distortion /JP'^ in both 1024^ and 
128^ simulations show critical exponents r and r'. Although apparently unrelated to the previous underlying critical 
exponent 77, this exponents r and r' quantify the fractality of the strain-history-dependent plastic distortion. Figure 18 
shows the correlation functions of the strain-history-dependent deformation "0, with the critical exponent r" close to 
2, which implies a smooth non-fractal field, shown in Fig. |16| ^c) and (d). All measured critical exponents are listed in 
Table mi 

Figure 17 shows the power- law dependence of the rotational invariants C^J^" and C^^ (they overlap). According to 
the definition (i^- = ik^ipj, we can write down the Fourier-transformed forms of Eq. (E2| and Eq. (E3) respectively 



2(/3S'''/3]^")(27r)3^(k) 



fcifcjV'j(k)'0i(-k), 
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V 



(E5) 
(E6) 



Except the zero-wavelength terms, the same functional forms shared by these two rotational scalars explain the 
observed overlapping power laws. 



2. Stress-stress correlation functions 



As the system relaxes to its final stress-free state, we can measure the fluctuations of the internal elastic stress 
fields, using a complete set of two rotational invariants of correlation functions 

CL(x) = «(xX(0)), (E7) 

Cj;(x) = (ajf (x)4f (0)); (E8) 

and in Fourier space 

Ct(k) = ^ajf (k)?jf (-k), (E9) 

C(k) = :^?if (k)aj'5*(-k). (ElO) 

Because is symmetric, these two correlation functions form a complete set of linear invariants under rotational 
transformations. 



3. Energy density spectrum 



The average internal elastic energy £ is written 



1 
V 
1 



2 iJ 



int int ^ ^int iiit 

ij ij 1 _|_ ^ " 33 



where, in an isotropic bulk medium, the elastic strain e° is expressed in terms of a" 



1 



^ , —int 



l + v 



We can rewrite Eq. (Ell I in Fourier space 

d'^k 1 



E 



1 



V J (27r)'*V 



5jf(k)5jf(-k) 



ajf(k)ai;;*(-k) 



Substituting Eq. (E9) and Eq. (ElO) into Eq. (E13) gives 



E = 



cfi'k 1 



cL(k)-— c?;(k) 
1 + 



(Ell) 



(E12) 



(E13) 



(E14) 
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lot (CGD),7 = -2.6S 
tr(CGD), Y = -2.65 
tot (GOD-MDP),7 = -1.65 
ti-(GOD-MDP),7 = -1.65 
tot (GOD-LVP),7 = -1.95 
tr(GOD-LVP),7 = -1.95 




5 



(b) 



le-6 




. tot(CGD),y = -3.1 
tr(CGD), 7 = -2.9 

. tot(GOD-MDP),y = 
tr(GOD-MDP), y=■ 



^tot(GOD-LVP),y= 
-- tr(GOD-LVP), y = 



0.02 



ks 



0.2 





O.lr 
(f) I 



«tot(CGD),y" = -l 
-- per(CGD),y" = -1 
^ tot (GOD-MDP), y" = -1 

► - per (GOD-MDP), y" = -1 
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FIG. 19: (Color online) Stress-stress correlation functions C"(k), elastic energy spectrum E{k), correlation func- 
tions of the stress- full part of GND density C'''^ (fc). Red, blue, and green lines indicate CGD, GOD-MDP, and GOD-LVP, 
respectively. All dashed lines show estimated power laws quoted in Table [Till 



If the stress-stress correlation functions are isotropic, we can integrate out the angle variable of Eq. (|E14|) 

M 



1 + 



CUk) 



(E15) 
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TABLE III: Power-laws relations among C^ik), E{k), and C (k). (d represents the dimension; P.Q. and S.T. represent 
'Physical Quantities' and 'Scaling Theory', respectively.) 

Simulations 

P-Q- S.T. Climb&Glide Ghde Only (MDP) Glide Only (LVP) 



20(1024^) 3D(128^) 2D(1024^) 3D(128^) 2D(1024^) 3D(128^) 

Ctot{k) 7 -2.65 -3.1 -1.65 -3.0 -1.95 -3.1 

Ctr{k) 7 -2.65 -2.9 -1.65 -3.0 -1.95 -2.9 

E{k) 7 + d-l -1.65 -1.1 -0.65 -1.0 -0.95 -1.1 

cC{k) 7 + 2 -0.65 -1.0 0.45 -1.0 -0.05 -1.0 

C^!^r{k) 7 + 2 -0.65 -1.0 0.45 -1.0 -0.05 -0.9 



where f{d) is a constant function over the dimension d 

fid) 



l/(87r) d^2, 
l/(87r2) d = 3. 



(E16) 



Writing the elastic energy density in terms of the energy density spectrum £{t) = E{k, t)dk implies 



l + v 



CUk) 



(E17) 



4. Correlation function of the stressful part of GND density 



According to Eq. (14), the stressful part of GND density is defined as 



Substituting Eq. (E12) into Eq. (E18) gives 



^ a / _int c _int 



The complete set of rotational invariants of the correlation function of includes three scalar forms 



Cf,"(x) 



{pt,{^)pfM), 
{pf,{^}pfm), 



(E18) 



(E19) 



(E20) 
(E21) 
(E22) 



where Cf^ (x) is always zero due to pf^ = 0. 
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Substituting Eq. (E19I into both Eqs. (E20| and (E21| and applying the Fourier transform gives 



CL(k) 



1 



2/i2(l + !y)2 \V 



4^2^ 



fc2 / I 



1 + 



(E23) 



(E24) 



where we ma ke us e of t he equihbrium condition diaij — and thus kiaij = 0. Substituting Eqs. (E9) and (ElO) 
into Eqs. (|E23[) and (|E24[) 



CL(k) 



_fc2_ 

4^2 



Ct(k) 



2t/ 



rC(k) 



(E25) 



^per (k) 



U2 r 



4^2 



CL(k) 



1+1/2 



rCr.(k) 



(E26) 



Here we can ignore the angle dependence if the stress-stress correlation functions are isotropic. 



5. Scaling relations 



According to Eq. (E17), the term k suggests that the power-law exponent relation between E and C is 

7' = 7 + d — 1. 



(E27) 



Again, both Eqs. (E23l and (E24) imply that the power-law exponent relation between C and is 

7" = 7 + 2, (E28) 

regardless of the dimension. 

Table III shows a nice agreement between predicted scaling and numerical measurements for power-law exponents 
of C"' , E, and C . These relations are valid in the presence of residual stress. 

During the relaxation processes, the elastic free energy follows a power-law decay in time asymptotically, seen in 
Fig. [6j All the above measured correlation functions of elastic quantities share the same power laws in Fourier space, 
albeit with decaying magnitudes in time. 
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Dislocations which cancel at the macroscale may be geometrically necessary at the mesoscale. See Sec. Ill for our rationale 
for not including the effects of SSDs (whose Burgers vectors cancel in the coarse-graining process). 

Changing the initial reference state through a curl-free plastic distortion (leaving behind no dislocations) will change /3'''^ 
but not P^'^; the former depends on the history of the material and not just the current state, motivating our nomenclature. 
For our simulations with external sheaPl, the k = of /Jf"'^ couples to the boundary condition. We determine the plastic 
evolution of the k = mode explicitly in that case. For correlation functions presented here, the k = mode is unimportant 
because we subtract P'^ fields at different sites before correlating. 

In real materials the dislocation dynamics is intermittent, as dislocations bow out or depin from junctions and disorder, 
and engage in complex dislocation avalanches. 

There are two uncontrolled approximations we make. Here we assume that the continuum, coarse-grained dislocation density 
p — determines the evolution; we ignore SSDs as unimportant on the sub-cellular length-scales of interest to us. Later, 
we shall further assume that the nine independent components of pij all are dragged by the stress with the same velocity. 
Only misorientation mediating dislocations are counted. 

More precisely, equation (4) of Ref. contains two different definitions for L''; the one in Eq. |51l and Lf^' = [{Vsig'^ij — 
ofij)^] = [i'^sQsfj'^)^]- L^' is of course a strain rate due to SSDs, but since varies in space U'' is not equal to . 
Ref. 1911 suggests using a two- variable version of the SSD density, p'^'^^(x, x') = — f>^(x), making the two definitions 

equivalent. 

In these analyses of TEM micrographs, the authors must use an artificial cut-off to facilitate the analysis. This arbitrary 
scale obscures the scale-free nature behind the emergent dislocation patterns. 



